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O ' Abstract 

o, 

CN , We test an optimised hopping parameter expansion on various Z2 lattice scalar 

P3 ' field models: the Ising model, a spin-one model and A(/>^. We do this by studying 

^ . the critical indices for a variety of optimisation criteria, in a range of dimensions 

and with various trial actions. We work up to seventh order, thus going well beyond 
previous studies. We demonstrate how to use numerical methods to generate the 
high order diagrams and their corresponding expressions. These are then used 
, to calculate results numerically and, in the case of the Ising model, we obtain 

pg I some analytic results. We highlight problems with several optimisation schemes 

^D ■ and show for the best scheme that the critical exponents are consistent with mean 

^ , field results to at least 8 significant figures. We conclude that in its present form, 

f~^ \ such optimised lattice expansions do not seem to be capturing the non-perturbative 

■^ ' infra-red physics near the critical points of scalar models. 

^' 1 Introduction 

Oh! 



00 



X 



The aim of this paper is to study the LDE (Linear Delta Expansion), a non-perturbative 
method. Finding alternatives to weak-coupling perturbation expansions is vital as there 
are many problems where such methods fail even when the theory is weakly interacting, 
e.g. at phase transitions. Perhaps the most accurate non-perturbative method is lattice 
5^ \ Monte Carlo, but this has trouble with fermions, non-zero charge densities and dynamical 

problems pOj. However, there are very few alternatives: large-N for continuum calculations 
is one example, a hopping parameter expansion on the lattice is another. Unfortunately, 
these non-perturbative methods can be poor e.g. see the expansions for critical exponents 
r\ and 7 in three dimensions using large-N approximations in j2|. Thus it is important to 
develop alternative non-perturbative methods. 

LDE is just such an alternative non-perturbative method. Basically, it is an optimi- 
sation of an existing standard expansion scheme. The earliest appearance of the term 
LDE we know is in jSlll] and some of this emerged from the logarithmic Delta Expansion 
work of jSl El 13 EI • However LDE is such a generic scheme that it has appeared in many 
contexts, has been rediscovered many times and under many names: Optimised Per- 
turbation Theory jHl E|, Gaussian Effective Potential approximation ^T], Variational 
Perturbation Theory [l2l UHl Cll CH], Order Dependent Mapping [2j, Screened Pertur- 
bation Theory ^HI, method of self-similar approximation jT2|, the variational cumulant 

*E-mail: T.EvansOimperial . ac .uk, WWW: http : //theory, ic . ac .uk/'^time/ 
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expansion ^H], action- variational approach ^H]- The method has been apphed to the 
evaluation of simple integrals [H EIH El 1221 12^] , solving non-linear differential equations 
[21 , quantum mechanics 011311111111122112312312111211121] and to quantum field theory 
both in the continuum P[I3[I21[I3I23I23E3EI1E21E3E3E3 and on a lattice 
(discussed below). ^ 

An important and recurring question is how good is LDE? It appears to be a quite 
arbitrary scheme though, as we will try to highlight in the next section, many of the 
supposed limitations of LDE are common to all other expansion methods. Likewise, 
widely accepted methods are not always as good as one might have supposed, see for 
example the large-N expansion of critical exponents in Pj. LDE has, in fact, been tested 
rigourously in QM where it has been shown to lead to fast convergence of calculations 
of "zero-dimensional" path integrals j23 E] and of the anharmonic oscillator ^l 123 
1^ . However, such rigour is probably unobtainable in QFT (see [23j for an interesting 
discussion), where in practice the only simple test is to compare results against physical 
results, or, failing that, Monte Carlo calculations. 

One of our goals was to look at the accuracy of LDE calculations on the lattice, 
where we are in effect optimising hopping parameter expansions. The method can go 
beyond Monte Carlo results, e.g. one can easily handle finite charge densities [36J, but 
here the idea is to focus on calculations where there are good Monte Carlo results, if 
not physical measurements, as these can then provide us with the 'true' answer against 
which to test our results. We chose to look at critical exponents in 0^ field theories 
as these are computationally accessible, they have been studied extensively theoretically 
and they have physical relevance, as shown by the available experimental data for critical 
exponents [T3] . 

Another goal was to extend the accuracy of existing LDE lattice studies. The liter- 
ature on lattice LDE contains work on scalar theories |HITH ll!?7|IHHl lHn i lin i liT | li ^ liH l 
113 113 Hi , pure gauge models |3II3l23l23E3ll3ll3ll3ll3and gauged Higgs models 
|47| HH] . Where they have made comparisons with other techniques, the conclusions are 
usually positive. For instance (Jones strong coupling in gauge theory) or the comparison 
of gauge Higgs models near transition [Tfl I48j . However, as far as we know, all the work 
in the literature, except one pure gauge study ^^, has worked only to third order in the 
expansion, doing the expansions by hand. To go to higher orders requires a computerised 
approach to the expansion. Thus our second aim is to investigate how this could be done. 

The codes we developed were easily adapted for studies of both the full 0"^ QFT and 
spin models. Thus we will also look at the Ising and spin-1 models which fall into the 
same universality class as the full 0^ model |3 El IHI- The feasibility of LDE for full 
field theories at higher orders will therefore be addressed while the numerical simplicity 
of the spin models can be exploited at various points. Specific results will be given for 
all three models. 

In the next section we will outline the LDE in the context of our 0^ field and spin 
theories on the lattice. Section three will examine our numerical results in detail, com- 
paring these with our more limited analytical results in section four. The final section we 
summarise our conclusions. 



^The wide applicability of this approach under different guises means our list in this paragraph is 
only a sample of LDE type papers, suitable as a starting point for further research. We do not claim to 
establish the precise history of the method. Also see P^ for comments about names. 



2 The Free Energy Density in LDE 

2.1 General Discussion of LDE 

We will start by developing the LDE in very general terms by considering a theory 
described by the full action 5* which depends on physical parameters p. We will be 
looking at scalar field theories where p will be masses and coupling constants. 

It is helpful to focus on the calculation of a specific quantity, though the principles 
are the same for all quantities, so consider the partition function Z 

Zip) := tr{e-^}. (1) 

The starting point for all expansions in QFT is that with the physical action for the 
theory of interest, S = S{p), there is no exact solution^ for Z. We therefore turn to some 
solvable theory described by an action 5*0 = So{v) which depends on a set of parameters 
V. These parameters v of the solvable action 5*0 will be specified later but for now we 
just note they may be different from the physical parameters p of the full action S. The 
idea is to expand about the trial action 5*0, exploiting its solubility. To do this we replace 
the action S by Ss the S-modified action, i.e. 

S{p) -^ Ssip,v) = Soiv)-5{Soiv)-Sip)) (2) 

Then rather than considering the complete partition function Z of (Q), we look at the 
6-modified partition function Zs 

Z{p) = tre~^ — y Zs{p,v) = tr{e~^^(^)} = tr{e~^»(^)e-'^('^«(^)~^(P))} (3) 

Clearly we need to set 5 = 1 at the end of the calculation to return to the appropriate 
action of the physical theory, lini^^i Zs = Z. Thus the 6 is merely a book keeping 
parameter introduced to keep track of the terms in the expansion about 5*0 



Zs = y2-rZ4p,v) (4a) 

^-^ nl 

n=0 

Z4p, v) = tr [e-^°(^) {So{v) - S{p)r] (4b) 



where we have defined Z„ as the n}'^ order term in the expansion of Zs. The expansion is 
in terms of increasing powers of 5 — 5*0 and the parameter 6 is just a parameter of technical 
use that records the order in the expansion of any given term in the calculation. Note 
that without further information, demanding |5| <^ 1 is neither necessary nor sufficient to 
control the convergence of this series, so we need not be concerned that 5 = 1 is required 
at the end of the calculation. 

Formally, summing all terms gives back the insoluble theory S when 6=1, but by 
assumption this can not be done here. However, we are free to choose 5*0 so that individual 
terms in the expansion, Z^, are obtainable, and in particular, the low order terms can be 
calculated in a reasonable amount of time. This is one major principle behind the choice 
of Sq. The idea is then to truncate the series, i.e. we consider 



^r=E> (5) 



71=0 



^In some rare cases in QFT, we may be trying to construct an exact solution via some expansion, e.g. 
as in two-dimensional sine-Gordon models. 



where the superscript of Z\ denotes the order at which the expansion is truncated. 

So far, all we have done is describe generic perturbation schemes. For later compari- 
son, let us consider standard (weak-coupling) perturbation theory in this context. Then 
S — Sq would be proportional to some small real coupling constant. For instance, in the 
0"^ theories studied here we would choose S — Sq = — J d'^xXcj)'^. Then we hope that for 
< A ^ 1 the small parameter A would ensure convergence, or at least that a truncated 
series gives good answers. In fact the former is clearly not true, QFT perturbation series 
do not usually converge, and the latter can also be false. The instability expected for A 
negative, however small, indicates that physical results will not be analytic in the com- 
plex A plane about A = 0. Thus a series in A can not be convergent. The behaviour of 
such series is hard to study in full QFT though it can still provide useful information by 
using ideas such as Borel summability.^ However, whether the first few terms are a good 
approximation, despite these issues, changes depending on the question being asked of 
the model. For instance, while weak-coupling perturbation theory gives accurate results 
for QED, at phase transitions weak-coupling perturbation theory is known to break down 
even in weakly-interacting theories. Ultimately then, all expansions in QFT are devel- 
oped using physical intuition that (S — Sq) is suitably small, confirmed only by comparing 
with physical results. It is only posthoc, and usually only in a limited way (if at all), 
that a rigourous mathematical basis for the expansion can be given. QED is a successful 
example of this process, and it shows that despite the mathematical uncertainties, the 
procedure outlined so far is worth following. So here, as in all expansions in QFT, a 
second principle guiding the choice of 5*0 should be that 5*0 is chosen so that terms up to 
n = R are sufficient to give a good approximation. Crudely speaking, {S — Sq) is to be 
thought of as "small" , while the size of the coefficient 6 is one and does not control the 
convergence. 

It is at this point that LDE departs from standard perturbative expansions, and 
where the non-perturbative aspects come in. In the LDE approach the parameters v of 
the trial action 5*0 are fixed using a variational method. These variational parameters 
V lie at the heart of the LDE method and provide the mechanism through which non- 
perturbative behaviour is introduced into the model. Before the expansion is truncated, 
upon substitution oi 6 = 1, the dependence on the variational parameters vanishes and we 
have the full theory which only depends on the physical parameters p in S, so Z = Z{p). 
However, a truncated series Z^ , will have residual dependence on these unphysical 
variational parameters v even when S = 1, i.e. Z^J^ = Z^J^(p,v). The final criteria for 
the choice of 5*0 is that it has some variational parameters v suitable for our subsequent 
exploitation. 

Clearly the dependence in the truncated answer on the arbitrary variational param- 
eters V is unphysical. Thus the final stage of the LDE method is to fix the variational 
parameters. Most importantly, this is to be redone every time we repeat a calculation at 
a higher order. Thus the values we assign to the variational parameters depends on the 
order at which we truncate the expansion so we are expanding about different theories 5*0 
at each order. Thus this procedure has been called an order dependent mapping 0. By 
contrast, if v was to be fixed at some specific order of the expansion and then used for all 
other orders, we would end up with nothing else than another perturbative expansion. 

The aim is to fix the variational parameters to values which will produce a result 

•^Rigourous results for LDE are possible in one space-time dimension when the problem becomes one 
of the Quantum Mechanics of an anharmonic oscillator. See also "53^ for a useful discussion of LDE and 
expansions in general. 



closest to the true physical one. The problem in achieving this goal is that there is no 
unique prescription which tells us how to do this. We know of five broad categories of 
methods used to fix the variational parameters. PMS (the principle of minimal sensitivity) 
and FAC (fastest apparent convergence) are the two most common. The names PMS and 
FAC, and their systematic definitions were introduced in [0]. However, prior to this we 
find a PMS-like criterion used in a study of a variational Hartree-type expansion applied 
to (j)^ field theory [211 and used in QM [23 12^1, and the FAC criterion used in jS2I 
to calculate the energy levels of an anharmonic oscillator by a variational perturbative 
expansion. In continuum calculations, the use of gap equations (choosing variational 
parameters such that the self energy is zero and the full propagator has a pole at the 
variational mass) is a third approach. A fourth is that of Meurice ^| who cut off the 
range of integration of the fields. 

However we shall focus largely on a fifth approach and fix our variational parameters 
by demanding that they minimise the free energy ^Sll33j. This choice is simply based 
on the physical principle that the free energy density f{p) of a system always tends 
to a minimum. Thus we seek a minimum of the function f^^\p,v) in the variational 
parameter phase space v whilst keeping the physical parameters p constant during the 
optimisation. This will fix the variational parameters to the optimum value, which we 
will denote with a bar as v — a function of the physical parameters p and the order R. 
Likewise the optimum free energy density will then be 

/(R)(p) ,= f(R)(p,v) < /(^)(p,«) Wv, (6) 

for a given set of physical parameters p. Assuming analyticity in v near v, we have 

Ifr - (7) 

At this optimal point v = v, 'small' variations in the components of v produce 'negligible' 
variations in /] , thus we are as close to being independent of the variational parameters 
as we will ever be. We will compare this minimum free energy principle against the other 
methods later in section 13.1.81 

2.2 Lattice scalar field models and LDE 

Now let us apply this to specific models. We will consider theories of a single scalar field 
on a hyper-cubic Euclidean space-time lattice"^ with lattice spacing set to be 1 and of 
various dimensions with action 

S = ^UL + 5'nul (8) 

^UL = 5^J0. + a0,2 + #f (9) 

ieA 

Snxjl = -'^X^ X^ ^i't^j (10) 

We have split the action into two types of term. S'ul contains the ultra-local terms, 
involving products of fields from the same lattice point. The non ultra-local terms S'nul 
have products of fields at different lattice points. The generic physical parameters p of 
the full action S are p = {k, J,a,g), though one can rescale the fields and remove any 

^See conclusions in section |5l for comments on hyper-tetrahedral lattices. 
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one of these parameters. A is the set of space-time points, (j)i is the field value at site i, 
and A/^"*" is the set of nearest neighbours to lattice point z in a positive direction. 

For a full field theory we choose 0j G M and exploit the freedom to rescale the field to 
set K = 1 so the physical parameters are then p f^ = (J, a, g). To make contact with the 
the usual notation used for the spin-1 model, we set (7 = 0, use the rescaling freedom to 
choose a = 1 so then Pspini = (^5 ^) ^'^d restrict 0j to 0.j G {+1, 0, —1}. For Ising models 
we further restrict 0j G {+1, —1} but we still have Pjsing = (fi^)- AH these models belong 
to the same universality class as the lattice 0"^ model [3 11311^1 • 

The first step in LDE is to chose the trial action 5*0 according to the three principles 
set out above, namely that 5*0 describes a theory that we can solve, 5*0 is a reasonable 
approximation to the full theory, and 5*0 contains non-physical, variational parameters. 

On a lattice, an action consisting of only ultra-local terms is soluble since the path 
integral factorises because the integral over the field at each space-time point is then 
independent of all the other field integrals. The insolubility of the full lattice 0^ action 
(jHl) comes from the non- ultra-local terms S'nul in the kinetic terms ^. Thus we must 
choose a pure ultra-local action for our trial action 5*0. Bearing in mind our second and 
third criteria for choosing a good 5*0, a suitable guess for 5*0 would be a trial action of 
the same form as the ultra-local part of the full action S but with arbitrary coefficients. 
Thus we introduce 

5o(«):=^Lo(0,,i;), (11) 

ieA 

where our trial action is ultra-local (solvable) and depends on the variational coefficients 
V. One choice is to make this of the same form as the ultra-local terms of the full 
Lagrangian, 0, to be a close approximation to the full theory. So we will work with 

Lo(0, v) := [j0 + k(P' + Icj)^] V := (j, k, I) (12) 

The variational parameters, j, k and I, are to be fixed by some optimisation condition 
once the expansion has been performed. In practice we do not always use the quartic 
variational parameter and in the Ising model both quadratic and quartic terms are trivial. 
In the end we will be left with expressions in terms of statistical averages {4'^}n. taken 
with respect to (fTTjl where 

(<^'>oW ■= ^ [d<P<P^exp{-Soiv)}=^f (13) 



Zoiv) := y rf0exp{-Lo(0,«)} = /o (14) 

Ip{v) := Jd^{(Prexp{Lo{(P,v)}. (15) 

Note that the ultra-local property of Lq is crucial and allows the factorisation into inte- 
grals over field values at a single site. 

This is a good place to note that the expansions we are interested in are traditional 
lattice expansions, and it is the variational aspect which is novel. If we were to fix our 
variational parameters equal to the appropriate physical ones, i.e. j = J, k = a and / = g 
in (fTTjl . then we would have an expansion in the the non-ultra-local term of 5* only. This 
is then a traditional strong coupling or hopping parameter expansion ^J. It is then not 
surprising that we can adapt the machinery developed for such traditional expansions 
[H EHl Ei E3 Eni EH EH for the general LDE problem. 



'S'^^ = 


- Si — (^i^NUL, 






^1 := 


= Sq + (52(S'uL ~ 


-So). 





2.2.1 Trick for ultra local terms in S — So 

If we recall the definition of the 5- modified action given in (j2)), we see that we have to 
expand exp{(5(5'o — S*)} to a given order and evaluate the path integral with respect to the 
ultra-local expl^o}. However, parts of {Sq — S) are also ultra-local and no more difficult 
to deal with than Sq. We can exploit this and make the diagrammatic evaluation below a 
great deal simpler. Thus let us make a further split of the action, replacing the physical 
action S by Sss 

S — ^ Sss = Si — SiSjsiUL, (16) 

(17) 

Li[0,] := {j + S,{J-j))^i + {k + 6,{a-kM + {l + S2{g-lM (18) 

The idea is that we can do an expansion in 6i first up to the desired order R, treating 
^i = 5*0 — 52('S'uL — 'S'o) as an intermediate variational action. The parameter 62 then 
counts the powers of — ^nul which contains all the non local terms in the calculation. As 
Si is ultra local, this expansion in 62 is just as feasible yet the diagrammatic methods 
used below are much simpler since there is only one type of term in — S'nul- This leaves 
us with an expression for a quantity of interest, such as Z, of the form 

n=0 

Zi4v,Si) = tT{e-'^'^P'^''^HSoiv)-Sip)r} (20) 

We can see that to return to the proper LDE expansion of ((21), we want to set 6 = 61 = 62 
in (fT?)|l and truncate the series rigourously to order R in 6. Thus we have to do a further 
expansion to order R — n in 61 oi each term ZiJyv). In general this would be doubling 
the work we have to do to produce an overall 8 expansion. However in this case we are 
expanding in powers of (S'o — Sul) about S'o and both are ultra local. This leads to this bi 
expansion being significantly simpler than the first expansion in 82 and the non ultra local 
terms. We will see this means that there is no need to use the full diagrammatic machinery 
a second time. This will be clearer when we have finished outlining the complicated 62 
expansion to which we now turn. 

2.2.2 Partition function expansion 

The expansions for all quantities can be expressed in terms of what we call statistical 
averages (Q),, given in terms of our intermediate variational action |T7j) : 

{Q)i ■■= ^y"p0Qe-^i (21) 

Zi := fv(t)e-^^ (22) 



For instance partition function can be written neatly as 

^^i = ^iE^((-^nul)">,. (23) 

n=0 

It is useful to introduce a diagrammatic notation at this point and to use the lan- 
guage of graph theory (e.g. see ^USH])- The whole space-time lattice is a graph whose 
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vertices are the space-time points. The edges of the lattice graph should be chosen to be 
those linking neighbouring space-time points which appear in non ultra-local products in 
the lattice action coming from the derivative terms. For our simple nearest neighbour 
approximation of the derivatives, the edges connect all nearest neighbour vertices. Thus 
each edge represents a unique nearest neighbour vertex pair and appears only once in the 
edge set of the lattice graph, i.e. our lattice graph is a simple graph. The non ultra local 
action can then be written as 

(24) 



S'nul - 


- >:>:■ 




i^^ j&Af+ 


i I - 


= ^(piCpj 



(25) 

From the diagrammatic point of view, performing the sum over i and j in S'nul is equiv- 
alent to running over all the elements of the edge set of the lattice graph. 

What this means is that the n-th order term in the 62 expansion of the partition 
function (j^^ is represented by a sum over all possible ways of choosing n edges from 
the edge set. We will call each choice a configuration. However note that we must 
include terms where we choose the same edge more than once. Thus in the language of 
graphs, our configurations include both simple and non-simple graphs. Only the simple 
configurations (i.e. ones which are simple graphs) are subgraphs of the lattice graph and 
these will be our main focus. We will see that the expressions represented by non-simple 
graphs represent straight forward generalisations of the expressions represented by some 
simple graph. We will use C to indicate a configuration of the lattice. 

Statistical averages have one important property given that we have chosen our vari- 
ational actions {Si here, 5*0 ultimately) to be ultra-local, namely the statistical averages 
of any polynomial of fields is also factorisable into terms depending on fields at only one 
lattice point. First the normalisation in our statistical averages factorises 



Zi = fv<j)e-^^ = H /"d0,exp{-Li(0i)} = {Z^f , 



(26) 



Zi := / d0exp{-Li(0)} (27) 

where Li is defined in (fTH|) . Note that for simplicity we are assuming here that the 
physical sources, p, are space-time constants, so that we obtain the same normalisation 
factor Zi whatever lattice site i is used in fl22|) (translational invariance). 
We can now write our statistical averages as^ 

{{<Pir{<P2r...{<p^r >i = ^y"i^0e-^H0ir(02)"^..(0.)"' (28) 

= {i<plr)^{i<p2r)l■■■{i<p^r), (29) 

(30) 



<-'ni <-'n2 '^Ui 



Jo Jo Jo 
and they depend only on the integrals over fields at a single space-time point 

Jniv,5i) := / rf0(0)-e-^i('^''^i'<^) (31) 



^The statistical averages {4>^)f. defined in Ijl^fl with respect to Lq of l|ll|) have identical factorisation 
properties but now in terms of the Ip integrals of (|15|l . 
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This gives us a simple translation from the statistical averages associated with configu- 
rations C to algebraic expressions: 

(^>i = «^" n t" ^32) 

t-t Overtax 

where Cvertex is the set of vertices in configuration C. ki is the connectivity of the i-th 
vertex, that is the number of edges in the configuration C which have one end at the i-th 
vertex. The configuration has n edges. 

It is important to note that since we have translational invariance, the expressions 
depend only on the connectivities of a given configuration and not on the precise position 
of the configuration on the lattice. It is therefore useful to consider graphs which are 
not embedded on the space-time lattice, and these we call diagrams, D. Every diagram 
is isomorphic to many configurations yet all these configurations represent the same 
algebraic expression. The number of different ways we can embed the diagram in the 
space-time lattice, that is the number of configurations isomorphic to a diagram is called 
the lattice constant of the diagram. For instance the expression for the partition function 
is then 

Zs. = Z,J2in5,r'-^cn ( J] ^) (33) 

Here the sum is over all possible graphs, simple and non-simple, connected and discon- 
nected, as these are the set of diagrams D. The co is the lattice constant for the graph 
D while n is the number of edges in the diagram. The last factor, bn, is a binomial factor 
coming from the fact that the same edge in a configuration of n edges, can come from 
any of the n factors of S'nul in the (— S'nul)" terms. Allowing for non simple diagrams, 
we see that 

e 

where the product is over all distinct edges, e, in the diagram, and n^ is the number of 
times each edge appears in the edge set of the diagram. Thus for any simple diagram, 6/3 
is always n\, while for a maximally non-simple diagram, i.e. diagram of n identical edges, 

2.2.3 Free energy expansion 

The partition function is not the easiest quantity to calculate. It involves all powers 
of the volume, here proportional to A^ the number of lattice sites. The free energy 
density, /, should be intensive and simpler to calculate. However, as the logarithm of the 
partition function it has a more complicated expression in terms of statistical averages 
so it is convenient to work with intermediate quantities called cumulant averages (also 
called cumulant expectation values or semi-invariants) and the expansion in terms of 
these averages is sometimes called a cumulant expansion or a cluster expansion^ . We can 



^The cumulant expansion owes its origins to studies in chemistry in the late 1930's. Starting in 1959, 
a series of papers [^13 ISl IE2| described the application of the expansion to the Ising model. Following 
papers jnHll^ further systematised the expansion, and applied it to the Ising and Heisenberg models. 
More recently, we find the cumulant expansion introduced and used in various ways and applied to a 
wide range of models. Some examples are j^l ^S BH ^| ■ The approach we take to presenting the 
cumulant expansion is motivated by j58j . 



define the cumulant averages to play the same role for the free energy as the statistical 
ones did for the partition function [53]. Adapting to our case, where we are keeping 
'^i('S'uL — 'S'o) in the exponential and expanding initially just in powers of — (52S'nul and 
truncating at order R we have by analogy with 



n=0 

/lo = -^In^i = -InZi, /i„ = - — ((-5'nul)")c (36) 

Note that we can use this to define the cumulant averages {| . . . ) . 

These cumulant averages have several important properties which we will use to cal- 
culate them. First they are linear so that we can use the diagrammatic notation 

( (-5nul)" >c = 5Z bDCD{D)^ (37) 

D 

where the sum is over all distinct order n diagrams, that is diagrams with n edges, and 
})£, and C£, are the binomial and lattice constant factors for diagram D exactly as we 
had above. Linearity together with the translation invariance in our model (p and v are 
space-time constants in our calculations) also means that every diagram is repeated once 
for each space-time point. Thus we can take advantage of this fact in future by calculating 
the embedding factors cd of the graphs per lattice site, at the same time discarding the 
factor of -^ sitting in front of the non-zero orders of the expansion in the above equation. 
However it is important to note that unlike lattice Monte Carlo calculations, our method 
works on a lattice of infinite size. 

A second property is the systematic link between cumulant and statistical averages 

(0e.---ee„>c= Y. (-l)'=-'(fc-l)!(ee,---ee„>^---(ee.---ee„>^ (38) 

partitions of (ei...en) ~^ 

K icictors 

where Gg must always be one of the operators whose sum makes Snul in (jTTjl . Thus 
here it is an operator </)j0j associated with one edge e = (i, j). It is vital to note that 
cumulant averages do not have a factorisation property analogous to that of the statistical 
expectation values p9|l . 

Finally, we come to the crucially important cluster property of the cumulant averages. 
Given the ultralocal 5*1 in the exponentials, the statistical averages in (jHH|l factorise if any 
of the operators Q^i does not depend on field values in any of the other Gg operators. One 
can quickly check that in this case the cumulant average vanishes. Put another way, the 
cumulant average of any disconnected graph is always equal to zero. A simple proof by 
induction can be found in [SEl, while the first general mathematical proof was presented 
in ISSI- 

This last property is of great importance to us, because it vastly simplifies our task 
of calculating the ( (5'nul)" )q terms. 

If we put all of this together we arrive at a final formula 



/^'^ := E;^/- (39) 



n=0 

/lo = -^In^i = - [In^i]^, /in = - E hr,'-^ m^]^_^ (40) 
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The sum is over the set V^ of connected diagrams of n-hnks. The notation [Q]m indicates 
that we have to expand the quantity Q to order m in a Si expansion 



iQL:^±' ''"^ 



,^0 ^'- ^(^^ 



(41) 

(5i=0 



A worked example is provided in appendix iBl 

2.2.4 Numerical evaluation of free energy 

The first step is to find the lattice constants cd- To do this we generate a complete set 
of simple connected configurations of up to R edges, modulo translation invariance. The 
algorithm used was inspired by one used in percolation theory [^11 IHZl and is outlined in 
appendix El 

We then have to see which of these configurations are isomorphic to a given diagram 
diagram D so that we can obtain the lattice embedding constant cd for each diagram. To 
do this we use a canonical labelling for our configurations, that is for each configuration 
we map it to an abstract graph that has a unique labelling. In this way we know that if 
and only if two configurations are assigned to the same abstract graph are they identical. 
To do this we used a small part of McKay's beautiful set of graph routines called nauty'^ 



It is a simple matter of combinatorics to generate the bf) and cd/N for all the non- 
simple but connected diagrams D up to order R by duplicating existing links in lower 
order simple diagrams.® 

Thus this stage of our computation provided a list of all connected diagrams up to 
order R, together with their bDCo/N coefficients for a specific lattice and in specific 
number of dimensions. However, the list can be used for any model where the only non- 
ultra local term is a nearest neighbour interaction. Our compiled programme producing 
this list of diagrams was relatively fast on, and used little memory of, a common desktop 
computer, typically taking less than a day. 

The next stages were implemented using the interpreted algebraic manipulation pro- 
gramme, MAPLE [nn]- This performed the transformation of the cumulant averages 
represented by the diagrams calculated previously (-D)„, into the expressions in terms of 
statistical averages, (...). It then had to expand each statistical average in terms of 6i. 
This involves replacing each Jp integral of ()31|) with the appropriate expression in terms 
of Ip integrals of (|T5|l using the Lq defined in (I12|l . Moving from Jp to Ip integrals requires 
explicit knowledge of the ultra-local terms in both the action and the variational action 
and so is highly model dependent. The complete expression for each diagram D is then 
truncated to order R in 6 = 6i =62. Putting the expressions for each diagram together 
with their coefficients gives an expression for the free energy. MAPLE then outputs op- 
timised C language routines for the free energy expressions. Thus each run providing a 
routine for a specific model in a specific number of dimensions and on a specific type of 
lattice. Results were checked by taking appropriate limits (easy in MAPLE) of code for 
the third order expressions produced independently in an earlier study [i5] . 

The advantages of using MAPLE were that development is relatively easy, it is ex- 
tremely easy to alter the code to suit different models or to produce code for derivatives 

^The name arises by forming an acronym from No AUTomorphisms, Yes?. 

^In our implementation, some of our non-simple diagrams were isomorphic to others in our list but 
this redundancy is relatively small. 
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of the free energy, and it produces fast C routines for the final stage of our analysis. 
The disadvantage is that as interpreted code it is relatively slow. Never the less, using 
only basic optimisations, the MAPLE code for seventh-order ran just within the memory 
(1Gb) and speed (one day) limitations of a desktop computer. It is clear that to go to 
higher orders a compiled version of these algorithms would have to be used. The main 
job of the MAPLE code is to produce and manipulate truncated polynomials with a rel- 
atively limited types of coefficient. Thus writing explicit compiled code for this part is 
feasible. 

The third and final stage of our numerical evaluation was implemented as another com- 
piled programme. Using the free energy routines provided by the MAPLE programme, 
the task was to find a set of values for the variational parameters v which minimised the 
free energy for a given set of physical parameters p. The Ip integrals required (for the 
full field theory case) are relatively straightforward to integrate, essentially dominated by 
a single peak in the integrand, so standard routines of [70] are sufficient. Our experience 
also showed that standard techniques to find the minimum in the multi-dimensional vari- 
ational parameter space [ZO], that defined by v, were successful. The practical problems 
encountered lay elsewhere. 

The first problem is that the expressions for the full seventh order free energy are 
enormous. The files for the free energy f^^^ grew as follows: R= 1 used 6Kb, i? = 3 used 
40Kb, order R = 5 used 0.5Gb and for i? = 7 we had a 5Gb file. We found we needed 
about 1.5Gb of memory to compile our seventh order routines. A rough guess for the size 
of the 9th order code file would be approximately 50 gigabytes in size, beyond the limits 
of desktop computing. 

However, this size indicates part of another serious problem, namely numerical accu- 
racy. Our 7th order code is a sum of a vast number of terms, each term being a product 
of several variational and physical parameters, some very large integer constants b£)Cd/N, 
and several Ip integrals. Each of the Ip integrals must in turn be evaluated numerically 
involving many additions. Worse near a critical point where typically j ^ there is a 
large variation in the size of the terms. Some, such as Ip for p odd, approaching zero 
as j -^ 0. Further, we are interested in variations of the free energy as its parameters 
are altered by small amounts. For instance, to locate the minima in variational param- 
eter space, we will often need to compare f^^\i = 0) against f^^^{j = e) for small e. 
Variations in the physical parameters are used to calculate physical quantities such as 
(^(py Thus tiny terms in the expressions for the free energy may well be significant for the 
overall change in the free energy. Just as bad, it is important to ensure that the variations 
in large terms, which may be only small fractions of the total value of that term, are also 
accurately calculated. A careful study of even the third order calculations suggests that 
the double precision accuracy available on typical C++ compilers for desktop PCs may 
not be sufficient. Thus we implemented our codes using arbitrary precision arithmetic, 
provided by the GNU multiple precision routines [71]. This enabled us to choose the 
accuracy of all our computations at the start of a run. For instance we were able to use 
256 bit arithmetic providing about 70 decimal places of accuracy. Further details and 
examples are provided in appendices O and O 

2.3 Physical Quantities 

The free energy is not the only quantity of interest. We can also calculate various deriva- 
tives of the free energy and it is the behaviour of these at the critical point which are 
directly linked to critical exponents. 
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Let us use the example of (</>) to illustrate issues which apply to all our calculations 
of derivatives of the free energy. Not only is this quantity the order parameter for the 
expected phase transition, but it is also the quantity which will give us an estimate of 
the critical exponent /3 (see equation ^7\i ). Analytically, (0) is related to the free energy 
through 

<^> - % (42) 

However, we are not dealing with the full free energy density, but rather with the trun- 
cated version /'-'^^ For this reason, we shall take a closer look at how we can make use 
of the above equation to calculate the field expectation value. 

One way of proceeding is to use the values of our optimised free energy f^^\p) := 
f^^\p,v) directly to numerically estimate <^0). After all, once we have optimised f^^\ 
we do have our best guess for the free energy and that is meant to contain all the 
thermodynamic information. The idea is straightforward to implement using 

This way, we have to make two optimisations of Z^^-*: once using some value J for the 
physical source, and another time using J + e, where e is some very small number. Note 
that this involves taking the difference of two free energies which may be very similar in 
size and we are only able to do this because we can choose whatever accuracy we require 
for our numerical calculations. 

The alternative approach is to produce explicit code for {^0) (p, v) and to set the 
variational parameters, v = v{p) to give a result 

The subscripts on the (j)^^^ indicate where the result was obtained via differentiation of 
optimised free energies (diff ) or via substitution of values of v which optimise the free free 
energy into the explicit expressions for (p^^^ (expl). This is straightforward in principle 
as it is relatively trivial to implement the derivatives required by P^ in MAPLE on the 
general algebraic expressions, and then it can produce optimised C routines for the final 
numerical evaluation, just as we did for the free energy. 

The difference is with 0jjg we must produce optimised variational parameters v(p), 
functions of the physical parameters p, at two physical values, J + e and J while for 0exp 
we only optimise the free energy at J (other physical parameters held constant). 

However, it is not quite so simple and it is instructive to have a more detailed look at 
this explicit code approach to the physical quantities. From equation P^ we see that 



^^^^ = 81 = -NZdJ ^^'^ 

where the generating function Z is of the form 

Z = tl e-'S'o(i,fc,0g'5A5(J,a,g;j,fc,Z) /^q\ 

with AS* = Sq — S. The second exponential in the above equation is the one that is 
expanded and truncated at some order R. Also in this term is the AS term that carries 
all the J dependence. Thus, by using equation (flKjl as a basis for calculating (0), i.e. 
differentiating / by the physical source J, we would bring down a power of 6 as well as a 
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factor of (j). If this was done on expressions for /^, then we would actually be producing 
an expression for (5(0) which has only {R— 1) terms in the delta expansion, and we have 
(0) only up to 0{5^~^). Put another way, the lowest order term in f'^^^ is In(Zo) which 
is independent of physical parameters like J and this the term lost. 

One solution is to develop a new diagrammatic expansion for (0) but this is costly 
and unnecessary. The solution to this loss of one order is to note that the first exponential 
contains 5*0 with a term linear in yet no factor of 5. So if we were to differentiate with 
respect to the variational source j instead of the physical source J, we bring down a sole 
factor of 0. The problem is, of course, that Sq in AS* carries j dependence too, so direct 
differentiation by the variational source would be of no use as it would produce other 
terms. Yet, all is not lost, because we can still use a trick to extract a factor of from 
5*0 without interacting with the AS" terms. 

We notice that the j coming from the e"'^'' term in the partition function is actually 
never seen explicitly in the full development of the free energy density — it is always 
hidden in the Ip integrals. On the other hand, the j coming from the e''^'^ term appear 
only in factors of 

AJ:=J-j (47) 

which are polynomial coefficients in the 61 expansion of the Jp integrals. It is easy then 
to keep 6 J constant while taking a partial differential with the remaining j factors by in 
effect only differentiating the Ip factors with respect to j 

Q- - -Vl ^4^) 

To perform the g. operation, we can go through the explicit expression of each /„ 
(c.f. equation (1121)) ^^^ replace every occurrence of Ip by the right hand side of equation 
(I48|) . It is straightforward to do the necessary replacements of the Ip in MAPLE and it 

produces optimised C code for (0) {p^v). Finally, for particular physical values p we 
can substitute the appropriate variational parameters t), which were determined by the 
relevant optimisation condition (minimising free energy, PMS, FAC etc.). Thus we can 

obtain (0)^^^^(p) = (0)^^^j(p,^). 

Another physical quantity of interest is the susceptibility x- It is defined by 

X-=-5^ (49) 

and is estimated in much the same way as (0) above. That is to say, we will use both 
the direct numerical approach, and the approach in which the individual orders of x^^'' 
are explicitly obtained. The direct numerical approach makes use of: 

AR) _ /(^)(J + 6)-2/(^)(J) + /(»)(J-e) 

Adiff — 75 l^uj 

Alternatively, finding the explicit cumulant expansion of x proceeds by applying the 
recipe we introduced for calculating (0): 

X-'«> ^ ^^^^^ (51) 



14 



with A J defined in fl47|) . Writing x*^^^ as an expansion, we identify the individual orders: 

xifpl(P;^) = X^-^Xn =^ Xn = ^^ (52) 

n=0 ■ -' 

As for the calculation of ((/)), we employ Maple to form optimised routines for Xcxpi(P5 '^)- 
The optimum value x^^^ is given by substituting the appropriate set of physical parame- 
ters and the values for the variational parameters obtained by optimising the free energy 

xSi(p):=xSi(p;^^). (53) 

To summarise, we have identified two different approaches for evaluating some phys- 
ical quantity Q (once the free energy density f^^^ has been used to fix the variational 
parameters). We can estimate Q directly using the optimised free energy density f^^\p), 
taking differentials of this function. Variational parameters lose significance once they 
have been used to optimise f^^^ (c.f. equation (021) )• Alternatively we produce an explicit 
expansion for Q^^\p;v). Some scheme is then needed to produce the optimised values 
of the variational parameters v, which could be derived either from the Q^^\p;v) ex- 
pression itself in which case the optimised values vary with the quantity Q of interest. 
Alternatively the values which optimise the equivalent free energy expression might be 
used. However the optimised values v are obtained, they are used to give the estimated 
result for Q, Q^^\p) := Q^^\p] v). One of our aims will be to produce results using both 
methods and to compare them. 

3 Results 

3.1 Ising Model 

We will start our discussion of results by looking at the Ising model since this model 
reproduces most of the behaviour spin-1 model and of the full 0^ model with the benefit 
of being numerically less complex. This numerical simplification comes from two aspects: 
the restricted degrees of freedom, 0j = ±1, and the reduced number of variational param- 
eters, V = [j). The trial action is defined with only one variational parameter because 
ultra-local quadratic and quartic terms are constant in the Ising model. The definition 
of 5* for the Ising model (jH)) leaves the model depending on two physical parameters 
p = {k, J). We will refer to k, as the inverse temperature^ . 

The restriction of the field to ±1 reduces the Ip factors to simple functions 



2cosh(j) p even 
-2sinh(j) p odd 



ip={ t:2::. r::: m 



3.1.1 Fixing the Variational Parameter 

Let us first consider the case with no physical source J. Then, as in |13], we find that 
the linear variational parameter, j, provides allows symmetry breaking to occur in this 
lattice LDE approximation. As long as J = j = 0, the trial action is invariant under 
the transformation of (pi —>■ —(pi. Once j ^ 0, the symmetry no longer holds, unless it 

^In the literature, the symbol /3 is predominantly used to denote the inverse temperature of an Ising 
model. We use k to avoid a notation clash with the critical exponent f3. 
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is accompanied by the transformation j —>■ —j. For instance the free energy for fixed 
physical parameters p = [k, J = 0) is an even function of j. Thus j = is guaranteed to 
be a turning point in variational space for J = problems. This behaviour is illustrated 
in figure ^ which shows plots of /^^•'(k, J = 0;j) vs. j for various values of k. The 
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Figure 1: Plot of f^'^^ (k, J = 0; j) vs. j for three different values. These are k = Kc — 10 
+ 10~^. The free energy is of the 3D Ising model at 7th order. 
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turning point at j = need not be the unique minimum but what we find is that the 
free energy profile has a simple form as a function of j. It is straightforward to apply 
our principle of lowest free energy to choose the optimal variational parameter value J. 
When K < Kc, the function /*^^^ indeed has a unique minimum at j = 0. For k -^ Kc, we 
find the region around j = flattening, and finally for k, > Hc the free energy displays 
two clear minima, either side of j = with the j = becoming a maximum. Both these 
minima are equally valid optimum points for j as guaranteed by the Z2 symmetry of the 
model ( J = 0). Numerically, the simplex optimisation method will choose one of the two 
minima, depending on the initial step of the routine. 

We can already make several deductions. The optimal value J is a continuous but not 
smooth function of k. This behaviour immediately tells us that our optimal value for the 
free energy will not be a smooth function of the physical variable n. Thus we clearly see 
that there will be some sort of phase transition at Kc and the variational source is itself 
a good order parameter even though it is not a physical observable. 

Also note that for the Ising model we find good global free energy minima for both 
odd and even orders. We will comment further on this later. 



3.1.2 Non-zero Physical Source 

For various calculations, we will require the behaviour of the optimised free energy as a 
function of the physical source J in the region of J = 0. Turning on the physical source, 
we break the Z2 symmetry and 'skew' the free energy vs. j curves. In the broken phase, 
one of the two minima acquires a lower free energy than the other, depending on the sign 
of the source. This effect can be seen in figure where we plot the free energy curves 
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with the system in the broken phase. The inverse temperature is k = Kc + 0.03, and 
we use three different values for the physical source. These are J = 0, J = 0.002 and 
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Figure 2: Plots of f^^\J]j) vs. j for a 3D Ising model, at 3rd order. The inverse 
temperature is k = Kc + 0.03, i.e. the system is in the broken symmetry phase. The three 
curves correspond to three different values of J . Each curve consists of 300 calculated 
points, which are connected by straight lines. The lines appear smooth due to the fine 
resolution of the data points. 

J = 0.004. With one local and one global minima, both in deep wells, we had to be 
careful to select the global minimum numerically. 



3.1.3 Critical Point 

The critical point, Kc, was identified by searching for a small interval where J(k_) = 
and J(k_|_) 7^ so that Kc = ^('^+ + ^-) i K'^+ ~ '^-)- Iii practice because the free energy 
is insensitive to changes in the variational parameter j (it is a minimum in j space) the 
minimum in variational parameter space is located to a much lower accuracy (in terms of 
the values for the variational parameters) as compared to the numerical accuracy of the 
routines. We need to know the variational parameters in order to find the critical point 
so uncertainties in these limit our knowledge of the critical point. This is discussed in 
detail in appendix iDl 

In table^we show the results obtained by searching for the critical inverse temperature 
at all orders of the 2D, 3D and 4D Ising models. Comparing these numerical values for 
the critical k, to those found exactly (see below in table El), we find the numerics were 
accurate to a fractional error of O(10~^^), i.e. to the accuracy quoted. This provides 
a further check of our numerics. A more detailed discussion of the errors in numerical 
identification of the critical point is given in appendix O 

In figure El we plot the values of Hc for all the orders. We also indicate the result 
Kc = 0.221654 obtained by Monte Carlo methods for the same model [Z21 IZHl HHj- The 
figure shows that if only the odd orders are considered, an imaginary line going through 
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Order 


Kc for 2D 


Kc for 3D 


Kc for 4D 


1 


0.2500000000000000 


0.16666666666666667 


0.1250000000000000 


2 


0.3333333333333333 


0.20000000000000000 


0.1428571428571429 


3 


0.3461538461538462 


0.20270270270270270 


0.1438356164383562 


4 


0.3768115942028986 


0.20963172804532578 


0.1464393179538616 


5 


0.3824833702882483 


0.21006903118305165 


0.1465228381635412 


6 


0.3944031482291211 


0.21343833354502731 


0.1475861410791982 


7 


0.3954564542557659 


0.21353497619553782 


0.1475854241672118 



Table 1: Numerical results for the critical points in the 2D, 3D and 4D Ising models, at 
all orders up to 7. Accuracy is at least to the last digit quoted (16th). 

the points seems to converge to a value close to the one predicted by Monte Carlo. The 
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Figure 3: The plot shows the values of k^ for the 3D Ising model, for all orders up to and 
including the 7th order. The Monte Carlo resuh is k^ = 0.221654 [7211731^9] . 

even points also show a smooth progression and appear to be giving reasonable behaviour. 
However, in slightly more complicated models the even orders do not give critical points 
at all. Experience tells us to distrust the even orders of our model and we will concentrate 
on the odd orders only. 



3.1.4 Field Expectation Value 

Consider first the evaluation of the expectation value via direct numerical differentiation 
of the optimised free energy, (0)jjg of (PSJ) . We present the results for orders 3, 5 and 7 in 
figure m We see that the expectation value displays similar behaviour for the three odd 
orders considered and qualitatively, the results are good. The expectation value clearly 
indicates the unbroken symmetry phase with a vanishing value. For n > Hc, we see the 
effect of spontaneous magnetisation. The only appreciable difference between the three 
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Figure 4: The plot shows {|0)^.jj(k, J = 0) vs. k ior R = 3, 5, 7, in the 3D Ising model. 
The physical source J = 0. The expectation value is calculated by direct numerical 
differentiation as given by equation ()43|). Each curve consists of 400 calculated points, 
which are connected by straight lines. The lines appear smooth due to the fine resolution 
of the data points. 



orders is the position of the critical point. Overall, the plot demonstrates that the LDE 
can be used to model spontaneous symmetry breaking with the expectation value as an 
order parameter. We shall look later at the near-critical behaviour, when we analyse the 
critical exponents. 

In figure El we have shown a plot of the expectation value calculated by optimising 
the free energy with respect to j first for a given k and J = 0, but then substituting 
the optimal variational parameter j into the explicit expansion of (^0) ,(/t, J = 0; j) to 



^(3) 



>(3) 



produce {(p)^^ j(/«, J = 0) (see equation (jSj)). In figure |31 we show a plot of {<p)^^ j and j 
for a range of k. These are results from a 3D, 3rd order run of the Ising model, again with 
the physical source J = 0. In broad terms, from the curve of (0) we clearly see again 
a phase transition at a point which we denote by k^, the critical inverse temperature. 
The critical points are the same for both approaches. 

By clearly acquiring a zero and non-zero value, depending on the phase, the expec- 
tation value is an order parameter of the model. Additionally, we find that the same is 
true for the variational parameter j. Thus we can use either order parameter to find the 
value of Kc by searching for the inverse temperature at which ^0) or j switch from a 
zero to a non-zero value. 

Note that the 3rd order curves using both methods look the same for both methods. 
In fact, accurate numerical comparison of the two methods reveal differences of ~ 10"^*^. 
Since we used e = 10~^° to calculate the numerical derivative (see equation pH|l ). this 
discrepancy is of purely numerical nature. ^° It is remarkable that the two methods 



""The same numerical calculation was done for the 5th and 7th orders too. 
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Figure 5: Plot of (0) and j vs. k for the 3D Ising model at 3rd order, with the 
physical source J = 0. Each curve consists of 300 calculated points, which are connected 
by straight lines. The lines appear smooth due to the fine resolution of the data points. 

produce the same results, especially considering how different the results will be for the 
susceptibility, discussed in the next section. 



3.1.5 Susceptibility 

As with the field expectation value, the susceptibility is another quantity that we can 
access by using optimised variational parameters or by direct numerical differentiation of 
the free energy density. The susceptibility is given by the second derivative: 



-(R) 

Xdis 



Q2jiR) fiR)^j + ^) _ 2/(^)(J) + /(«)(J - e) 



9J2 



(55) 



where the right hand side defines the direct numerical differentiation approach. The other 
approach is based on a function xLpil'^' J'-ij)^ which will assume the 'correct' value once 
an optimised j is substituted, i.e. x^^\n, J) = x^^Ki^^ J]j)- 

Figure El shows a plot of four different x^^\n, J = 0) curves. Three of those are 
Xexpi('^? J) calculated by optimising j first from the free energy, and then substituting 
this optimised value into x^^\ ^oi R = 3, R = 5 and R = 7. The fourth curve shows 
X^jg(fi;, J = 0) as calculated through direct differentiation (j3^ . This plot leads to impor- 
tant conclusions about the methods used to calculate each curve, and we shall discuss 
these in what follows. 

The curve denoting the direct numerical differentiation approach Xdis encapsulates 
the correct qualitative behaviour of the susceptibility. It is the case that x blows up at 
the critical point, due to correlations extending over the whole lattice (system). This 
is due to the fact that this approach is capturing the non-smooth nature of the field 
expectation value apparent in figure (jlj). The xLpil^' '^) form has only finite peaks at all 
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Figure 6: The plot shows xlxpii'^y J = 0) vs. k for i? = 3, 5 and 7 in a 3D Ising model. 
These three curves are calculated using a previously optimised j. The 'i? = 3 (diff)' curve 
is Xdiff('^) J = 0) calculated using direct differentiation, as given by equation (j33j) . Each 
curve consists of 200 calculated points, which are connected by straight lines. The lines 
appear smooth due to the fine resolution of the data. 



three orders, with the peaks getting larger as the order increases. The fact that the peaks 
shift towards higher k, is purely due to the value of Kc shifting at each respective order 
(see table Hj). It is easy to see (mathematically) why the peaks in the forms Xexpi('^' ^) 
become more pronounced at higher orders. The explicit expressions for x^^K^^ J'^J) ^"^^ 
a sum of a finite number of terms, each a product of Ip factors, physical and variational 
parameters. Like the expressions obtained in a high temperature expansion on the lattice, 
these expressions can not become infinite. Physically, this is a reflection of the fact that 
at, for example, 3rd order, our method can only correlate lattice sites spaced no more 
than three sites apart. Diagrammatically, at 3rd order, two sites i and j, sitting three 
sites apart from each other can be correlated only through the diagram: 



(56) 



If i and j happen to be further apart, the method will never be able to correlate them. 
The situation for R = 5 and i? = 7 is analogous: the method can only correlate lattice 
sites which are within R sites of each other and they can never see the infinite correlation 
length near the critical point. This is the reason why the xlxpi peaks grow as orders 
increase, but still remain finite for finite R. 



3.1.6 Critical Exponents 

The f3 critical exponent is defined in terms of the reduced inverse temperature k^ 
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Let us take the three-dimensional order seven results as a prototypical example. Using 
data produced with 155 digit accuracy, we fit 
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(58) 
The best fit 



using 30 digit accuracy on the fitting for 100 values of k G {nc, i^c + 10 ^°] 

is for /5 = 0.5 + 3 X lO^^^ and {n^Kcc^^t - 1) = 1-4 x lO^^o with x^ = 8 x lO-^o as figure 

m shows. The reference value Kc,exact is given later in table El The fit is very good so 
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Figure 7: Plot of the standard statistical measure x^ for the best fit of (jHH|l to 
the behaviour of the 3D Ising model 7th order results near critical point for a fixed 

(Kc/k, 



c, exact 



the fact that /3 is not exactly a half looks significant at first given the high calculational 
accuracy. 

However the fractional differences between this best fit curve and the data points 
are shown in figure |H1 This shows a small but significant systematic error that we have 
not been able to eliminate, perhaps due to a poor choice for the form of the fitting 
function (J3H)) . There is no reason why our approximate form should have exactly the same 
functional form as the full solution and we are working here to high levels of numerical 
accuracy. Fitting data in the same way but in the ranges k G {kc, Kc + lO"*"] for r = 6 and 
8 reveals that the all the errors and the deviation of the best estimate for (/? — |) reduce 
by two orders of magnitude as we use data ranges with r = 6, 8 and then 10, closer and 
closer to the critical point. If we fit a pure power law (c = in (J3H|) ) then the fit is slightly 
worse, the beta is different and again differs by the same order of magnitude from 1/2. 
Attempts to fit forms for other corrections to (J57|l to the data produce similar results 
but no consistency in any shift of /3 from one half. Thus, despite our high numerical 
accuracy, we conclude that the deviation from /5 = | is not significant and that our data 
is consistent with /3 = 1/2, again to about nine significant digits. 

Power law behaviour of the susceptibility near the critical point defines the 7 critical 
exponent as 

X OC \Kr\ ^ (59) 

To find the value of 7, we used the susceptibility Xdis calculated by direct differentiation 



of the free energy density (j33|l . 



Following our experience with (3 we merely fitted the 
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Figure 8: For the best fit of tlie critical behaviour of the 3D Ising model 7th order 



results (J = 0), a plot of ln((0)^ V((0> ) vs. ln|fi;^| for k G (k^, k, + 10"^°] with 
/3 = 0.5 + 3 X 10-10 and {KjKce^^ct - 1 
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simple power law form ()59|) to the data for In 

K G (kc, k,c + 10-*^]. The resuh is 7 = 1 * (1 ±10 , 

Finally the 6 critical exponent is defined by the following power law 



with K, taken in the broken phase and 
^) for orders 3, 5 and 7. 



(0>~|JP 



for Kj. = 



(60) 



Note that the S critical exponent has nothing to do with the bookkeeping parameter 6 
defined in the LDE methodology. Again we fitted the power law form (jUUj) to the data 
for (0)^jjj with K = Kc- The result of the fit at orders 3, 5 and 7 is 5 = 3 * (1 ± 10~^), to 
similarly high accuracies as for the other critical exponents. 

A summary of the critical exponents for the Ising model, as found by our LDE method 
is shown in table El For comparison, we give the values of the same critical exponents 



Exponent 



7 
6 

S/il/P + l] 



LDE 



0.5 
1 
3 

1 



MC+HT 



0.3485 
1.3177 
4.780 
1.000 



Table 2: The values of the three critical exponents (3, 7 and 6 for the Ising model in 
three dimensions. First as obtained by our LDE study (at orders 3, 5 or 7) secondly by 
a combination of Monte Carlo and high-temperature expansions |57| 1^ . 



as found by [58], using a combination of Monte Carlo simulations based on finite-size 
scaling methods, and high-temperature expansions. We note that the values obtained by 
the LDE match those predicted by mean field theory up to errors. 

One final test is to check the universal relationship 6 = 1 + (7//9) but as mean field 
values satisfy this, our values also satisfy this requirement within the accuracy of our 
measurements. 
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3.1.7 The Ising model in other dimensions 

So far we have given a detailed description of our results for the 3D Ising model. In table 
El we list critical inverse temperatures for all orders of the 2D, 3D and 4D Ising models. 
Again these are exactly solvable so can be used to check our numerical calculations. The 
critical exponents can then be measured as already described. Again they are identical 
to the 3D results, that is mean field values at orders 3, 5, 7 within the accuracy of our 
calculations. This independence of critical exponents with dimension is also characteristic 
of mean field theory but is not seen in the true results. 

3.1.8 The Ising model and other optimisation schemes 

Once the complete expressions have been obtained for the various scalar field models, it 
is relatively quick to compute them for a range of values in the Ising model with its single 
variational parameter and with Ip 'integrals' being given by simple functions. It is an 
ideal place to illustrate other aspects of the LDE method. We will just study the J = 
case. 

First let us look at optimisation scheme known as FAC - fastest apparent convergence. 
In this scheme one chooses the variational parameters such that the last term in the delta 
expansion (or sometimes the last r terms) is zero. The idea is that in well behaved series 
the last term can often be used to get an idea of the error due to truncation of the series 
so the 'optimal' result is where this is zero. 

The modulus of the fractional contribution to the free energy from each order in the 
full delta expansion is plotted in figure for the seventh order calculation. The spikes 




Figure 9: The eight terms in the delta expansion of the free energy of the Ising model 
in three dimensions up to order seven with J = and k = 0.23061777. Normalised with 
respect to the total free energy, the log of the absolute value is plotted. The numbers 
on one end of every curve indicate the order of the term plotted. Note that most spikes 
represent points where a term is going to zero but numerical and plotting limitations 
mean the curves are of finite extent. Likewise, the curves ought to be exactly symmetric 
about j = but near the spikes numerical limitations prevent this appearing on the plot. 
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ought to be going to minus infinity on this log plot as these are points where a term goes 
through zero. 

It is clear that the odd terms always go to zero at j = 0. This should be expected 
as all the odd integrals Jp are zero when j = 0, and at odd orders, for a trial action 5*0 
made up of only odd powers of fields, every contribution in the sum of such contributions 
making up the odd order expression must contain at least one Ip with p odd. Thus only 
an FAC procedure based on odd terms alone, will provide a suitable unbroken j = 
solution. Thus let us further confine ourselves to an FAC procedure where we demand 
that the optimal variational parameters v are such that the last term in the expansion is 
zero 

/Ffc(P):=E^/n(P,^), fR{p,v)=0 (61) 

n=0 

Now we notice that, at least for the value of k, used in figure Q, there are in fact 
non-zero values of j which satisfy the FAC optimisation criterion (jfilj) . There are four: 
two positive j and the same solutions with opposite sign (as there must be under the Z2 
symmetry at J = 0). Further investigation shows that these two new solutions appear 
first at KpAcc ~ ^c * (1-048 ± 0.002) where k\j ~ 0.2135 is the critical k value found 
in the same model but using the minimum free energy criterion ((7j). Now the problem 
for K > K,p^(^ fj is which minimum to choose — we have three distinct possibilities after 
symmetry is taken into account. This illustrates a general problem when using criteria 
such as FAC and PMS. They often present multiple possible solutions for the optimal 
variational parameters and one must use further criteria, sometimes no more than physical 
intuition, to choose one solution. The minimum free energy method we use guarantees a 
unique answer, up to symmetry, except at transition points. 

To finish with FAC, let us choose the j solution that has the lowest free energy as 
well as satisfying ()6H) for this seventh order J = example. As figure (jH)) shows this is 
the smallest non-zero j solution (the free energies are negative, so this has the largest 
free energy ratio when compared to the j = solution). This solution appears at about 
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Figure 10: The optimal non-zero j solutions (left hand plot) and the resulting free energy 
normalised to the free energy at j = (right hand plot) for the Ising model in three 
dimensions at order seven with J = 0. Note the smaller j solution (crosses) has the most 
negative free energy. 

J ~ 0.85 at K = Kp^(^ (J and so this FAC criterion is giving us a first order transition. Any 
physical quantity we calculate will suddenly change in value as we increase k through 
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(7) 

'^FAcc ^^ ^^^ i parameter is suddenly changed from to around 0.85. In some models 
this is a good thing. For instance in the Electroweak model for some parameter ranges we 
expect a first order transition and this sudden change in the optimal variational parameter 
solution can give this behaviour in LDE as obtained in j?7l|lH]. However, we are expecting 
a second order transition in this model. Our conclusion is that FAC is an unsatisfactory 
optimisation scheme for the optimised hopping parameter expansion of the Ising model. 
We next turn to PMS. In this scheme, one looks for turning points in the quantity of 
interest. Thus for a fixed set of physical parameters p the optimal values chosen for the 
variational parameters v will vary as we study different quantities. This is a standard 
procedure when PMS is applied to issues such as scheme dependence in particle physics 

[min]. 

First let us apply PMS to the free energy. The plot in figure Q shows how this works. 
For K < Kc there is a single turning point at j = at the global minimum of the free 
energy so the PMS gives the same answer as free energy minimisation in the unbroken 
phase. However for k > k^ there are three turning points and in principle PMS does not 
distinguish between turning points in the variational parameter space. In practice one 
often chooses the flattest turning point (for least sensitivity to unphysical parameters) 
but here there is little to choose between them. Failing this, with PMS one often falls 
back on choosing the turning point which makes the most physical sense. Here that might 
be the ones with the lowest free energy or with j ^ since we are expecting symmetry 
breaking. Of course with such additional arguments one selects the same solution as we 
had when simply searching for the minimum of the free energy. In this sense PMS is 
working here as well as the minimum free energy criterion. What it does highlight is that 
PMS need not give a unique solution nor a single solution. 

More interestingly we can try the PMS on other quantities. Let us look at the 
expectation value of as a function of variational v = {j} and physical parameters 
p = {k, J = 0}. This is shown in figure ITTl The main point to note is that (0) is an odd 




Figure 11: ((/))(«;, j) in the three dimensional Ising model at order seven with J = and 
fixed K. The labels indicate that {k - k^) = -0.05 (1), 0.0 (2), +0.03 (3), +0.05 (4) and 
+0.1 (5). Note that the new turning points appear at a kappa about 14% higher that 
Kc, the value shown with curve three. The curves are all odd in j so negative j are not 
displayed. 
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function of j and for large ranges of k, the low kappa regions, there are no turning points 
in (j). For larger (^'s there are four turning points, with one nice shallow local minimum 
as a function of j (maximum for negative j) appearing at some k ^ Kc- However it is 
clear that the PMS fails to give reasonable qualitative behaviour for this quantity. 

Thus for the Ising model, the FAC fails to work at all, while the PMS fails when applied 
to some quantities, and works on others^^ if implemented with some extra reasonable but 
ad-hoc criteria. On the other hand the unambiguous minimum free energy principle 
gives great quahtative behaviour for all quantities (when calculated using derivative of 
the optimised free energy) but it produces poor critical behaviour of a type seen in the 
mean field approximation. 

3.2 Spin-1 Model 

The spin-1 model differs from the Ising model in that can acquire the values ±1 and 
0. This means that a quadratic ultra-local term is no longer redundant (as it was for the 
Ising model), so S and 5*0 are (c.f. equation (jH))): 

So = 5^[j0. + A;02] (63) 

ieA 

Thus the Jp integrals ()15|1 are again simple functions 

r 1 + 2e-^ cosh(j) (p = 0) 
-^pO) ^) •= \ 2e~'^ cosh(j) [p even and positive) (64) 

y —le'^ sinh(j) [p odd) 

However, now we are optimising in the two-dimensional variational space oi v = {j,k). 
We set the parameter a = In 2 ^[ EO] . 

The optimisation with respect to two parameters causes no additional problems at 
odd orders. The j parameter acts as an order parameter as before. At low values of k 
k = 1 and j = gives the lowest free energy, i.e. where the variational parameters equal 
the physical parameters v = p and the trial action exactly equals the ultra-local part of 
the physical action. However, the gradient in the free energy with respect to j slowly 
decreases and becomes zero at some point and this is the critical point. For higher n 
values there are a pair of minima at j 7^ and k ^ 1. The free energy as a function of 
the variational parameters is shown for the third order in figure (fT^ for values of kappa 
either side of this critical point. 

Table 01 summarises the critical inverse temperatures found for the 2D, 3D and 4D 
cases. Exact expressions can be obtained but they involve integer powers of 1/e and are 
not given here.^^ 

The critical exponents behave in the same way as for the Ising model: /3 = 0.5, 7 = 1 
and (5 = 3 to over eight significant digits, irrespective of order or dimension. 

We present only odd orders, as the even orders do not produce good behaviour in 
the variational space. It is the behaviour in the quadratic variational coefficient, k which 
leads to the even orders failing to produce good minimum energy values. For instance 

^^ Presumably PMS works for those quantities even in j. 

^^This is because in the unbroken phase, the coefficient of the quadratic term k prefers to be equal to 
that in the physical action, here chosen to be 1. Thus the Ip integrals contain 1/e factors. 
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Figure 12: Contours of constant free energy (highest values at edges of plot) for the spin 
one model in three dimensions at order 3 for n = 0.9k,c (left) and k, = l.l/tc (right). Here 
Kc is a genuine critical point. 



Order 


K, for 2D 


Kc for 3D 


Kc for 4D 


1 
3 

5 

7 


0.5000000000000000 

0.5753424657534247 
0.6060227414719480 
0.6475789767717944 


0.33333333333333333 

0.36464088397790055 
0.37172944976396613 
0.37519182797115188 


0.2500000000000000 
0.2670623145400593 
0.2697296215999957 
0.2700940775993100 



Table 3: List of all critical points for the 2D, 3D and 4D spin-1 models, at all odd orders 
up to 7. The values are obtained by using the linear and quadratic variational parameters, 
j and k. The result obtained by Monte Carlo methods for the 3D case is Hc = 0.383245 



in figure (fT^ the free energy of the spin one model at second order is shown. The free 
energy behaves reasonably as a function of k and we can identify a value Kc where where 
the second derivative with respect to j is zero at j = 0,k = 1. However, there, as for 
all values shown the free energy has no turning point with respect to k and no obvious 
minimum. 

Interestingly, we tried reducing the number of variational parameters to one, namely, 
we used j only. The results produced in this way were the same as with the k variational 
parameter included, for all cases. The only differences occurred at even orders, at which 
we could not find any reliable results anyway. 



3.3 (f)^ Model 

As a reminder, the action and the trial action of the 0^ model are given by 

S = - ^ ^ </.,</,^. + ^ [ J0, + a0^ + ^04] 

Lo{(pi,v) := [j^i + k(t)^,+l(t)t] 



(65) 



(66) 



Thus the Ip integrals (fTH|) are now non-trivial and we are optimising in the two-dimensional 
variational space oi v = (j, k) if we set I = g or more generally we work with a three 
dimensional variational space v = (j, k, I). Note that here the inverse temperature k has 
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Figure 13: The free energy of the spin one model in three dimensions at order 2 at k/k^ = 
0.9(top blue curve), 1.0 (middle green continuous curve), 1.1 (bottom red dashed curve). 
At order 2, k^ is the point where the second derivative of the free energy with respect 
to i is equal to zero at j = 0, fc = 1, as the left hand graph indicates (its for k = 1). 
The figure on the right is at j = and shows that it is the lack of a turning point in 
the k variational parameter which prevents the method working at order 2. (Tim Maple 
analysis) 

been scaled to unity when discretizing the 0^ action (c.f. section l2^ . Instead, we use 
a, the lattice mass parameter, as the physical basis for 'driving' the system through the 
phase transition, so we have physical parameters p = {J,a,g). We set g = 25, to agree 

with ngcHi. 

The crucial difference in this model as compared to the Ising and spin-1 models, is 
the number of degrees of freedom of the field. For the 0^ model, the field is a continuous 
parameter, and the trace in the partition function becomes an integral. Therefore, the 
numerical calculation of the Ip factors is a computationally intensive task as they are 
now integrals whereas they were standard functions for the previous two models. We had 
to make sure that the accuracy of our numerical integration was appropriate given the 
accuracy problems encountered in manipulating the long expressions the LDE gives for 
the quantities of interest. For instance, when 256-bit arithmetic was used (which gives 
approximately 77 decimal places of accuracy) , we ensured that the integrals were always 
evaluated to a precision of at least 70 decimal places. 

Table 0] summarises the critical values of a at odd orders for the 2D, 3D and 4D cases. 
The critical values were calculated using only the linear and quadratic variational param- 



Order 


ac for 2D 


ac for 3D 


ac for 4D 


1 
3 

5 

7 


-14.585428484653014 
-18.407621223877752 
-20.106877371919612 
-21.240866931195218 


-10.007548118296870 
-11.560051936952422 
-11.918151992492446 
-12.138250831659820 


-6.956574606584720 
-7.887721708383015 
-8.018314178185453 
-8.062177866626536 



Table 4: List of critical ac for the 2D, 3D and 4D 0^ model, at all odd orders up to 7. 
The values are obtained by using the linear and quadratic variational parameters, j and 
k. 

eters. The critical exponents P, 7 and 6 turned out to be 0.5, 1 and 3, respectively, again 
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Figure 14: The free energy of the spin one model in three dimensions at order 2 at 
K = 0.9kc and k, = 1.1 Kc- Here Kc is the point where the second derivative of the free 
energy with respect to j is zero at j = 0, fc = 1. It shows that it is the behaviour with 
respect to the k variational parameter which stops the method working at order 2. (Tim 
Maple analysis) 

reproducing the values found in the Ising and spin-1 models to at least eight significant 
figures. 

The 3rd order result in 4D successfully reproduces and refines the result ac = —7.88 
found by a similar variational cumulant approach in [18^ . Two other studies of the 4D 
theory, based on Monte Carlo methods [ZSl IZEI found ac = —8.2 and ac = —8.275, 
respectively. A crude extrapolation of the trend of our 4D results, as given in table El 
gives a critical value of ac ~ —8.08. 

Remarkably, even in this model, the results using the linear variational parameter 
only are no different. The critical mass parameters ac, at all orders and in all dimensions 
are the same. We also found that including the quartic parameter made no difference to 
the critical behaviour. 

4 Exact Analysis of Ising Model case 

In this context, exact means no numerical evaluations were used. The approximate 
algebraic forms produced for the free energy were solved algebraically, including the op- 
timisation aspect. This is possible only for Ising model where there is a single variational 
parameter j. 



4.1 First order 

The lowest order, R = 1, LDE approximation to the Ising model is easily calculated. 
The derivation also provides a simple concrete example of the LDE method. For the free 
energy we find that 



-\nZo-{AS)^ (67) 

= -\nZo-Kd{<P)l-{j-J){<P)^ (68) 

We are using statistical averages with respect to the Sq action, defined in (fT^ so that 

Zo = Io, (0^>n = T^ (69) 
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They also factorise in the same way as the (Q), expectation values ()29j] . Note the absence 
of quadratic and quartic terms from (^AS')„ because of the simple form of the action and 
trial action for the Ising model. To find the minimum, we need to differentiate with 
respect to the sole variational parameter j, and set the resulting equation to zero. Using 
dlHl) we find that 

?fl^(2.i(4,\ + j-j){(4'\-{4')l) (70) 

We set the above equation to equal zero, and provided (0^)q 7^ ('^)o' ^^^ optimum j 
value, J, is found by solving 

j = J-2Kd{<P)^ (71) 

The equation is generally transcendental since (0)„ depends on j. 

For the Ising model, using the form of the Ip given in ()54|1 . it is possible to write this 
equation as 

j = J + 2Kdtanh{j) (72) 

In general this equation has three real solutions. If we wish to use the principle of 
minimum free energy we have to pick a solution which actually produces a minimum. 

However, one can quickly see that for J = then there is only one solution in the 
region k < Kc = \/{2d) and that is j = 0. Further, d'^f^^ydj'^ = at k = Kc where the 
j = solution turns from a global minimum into a local maximum (as happens at higher 
odd orders, see figure H]), and the correct value of J is a non-zero one for k, > Kc- 

One can now see what the critical exponents are. First let us note that the delta 
expansion for (0) in the Ising model gives to lowest order 

W = Wo + 0--^)((<^'>o-(Wo)')+^-((<^'>oWo-(Won (73) 

= - tanh(j) + ((j - J)+dK tanh(j)) (1 - tanh^(j)) (74) 

Near the critical point we expect j{J) to be small as J — ^ 0, k — + k^ = l/{2d). Normally 
we consider the optimal variational parameter j to be a function of the two physical 
variables k and J . To solve analytically, however, its easiest if we express the relationship 
the other way round, i.e. making expansions in terms of j and express the one physical 
parameter we vary in any critical exponent expression in terms of j. Thus keeping k = 
Kc and studying J{j), f{J{j)), (0)(j) in terms of a small J expansion, we find that 
J = J'^/3 + . . . and (0) = — j/2 + . . . showing that the critical exponent S of (jHIH) is 3. 
Likewise, for J = we find that {k/kc — 1) = J^/3 + . . ., (0) = — j/2 + . . . and so the 
exponent /? of (j^ is 1/2. 

The exact same exponents are obtained in an analytical mean field analysis. So is 
the LDE using minimum free energy optimisation the same as Mean Field theory? The 
answer is not exactly. In mean field theory for the Ising model on a simple hypercubic 
lattice we would write 

Z = /"D0exp{«:5^ 5^ 0,0,-J5^0,} (75) 



^MF 



/ Dt] exp{KNdv'^ - JNv + 2dvK ^Vi- J^Vi} (76) 

•^ igA ieA 

2Zocosh{2dKV-J) (77) 
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where rji = (f)i — v and v is defined, self-consistently, to be the expectation value of the 
field, i.e. 

v= ^ \,T = ta.nh{2dKv - J) (78) 

Comparing the variational parameter here, v with j of the our LDE calculation using 
lowest free energy optimisation, (f7^ . we see that the two are related by 

3 = J - 2dKv (79) 

However, the forms for the expectation values (0) in the two calculations are clearly 
different as comparing ()74|1 and ()78j) . using ()79j) clearly shows. The difference though 
becomes negligible near the critical point which explains why the two approaches give 
the same answer for critical behaviour. 

There is another way to view the link to mean field and that is to realise that mean 
field is merely another optimisation scheme within the LDE family. If we chose our 5*0 
to be 

5o = 5^Lo(0i,t;), Lo(0,t;)=fi+j0 (80) 

rather than that of ()12j) we see that the LDE parameter j is playing the role of a mean 
field, and comparing (jHOI) with (fTTjl . we obtain the relationship (f7n|) but for all j values, 
not just the optimal values. If we used an optimisation scheme where 

j:=2dK{(p)-J (81) 

then we would obtain MP as the zero-th order in this LDE scheme. 



4.2 Exact results for k, 



c 



The Ising model can be solved exactly within the LDE approximation with minimum 
free energy optimisation at least at the orders and for the dimensions studied here. This 
was done using the algebraic forms for the diagrams produced by our nauty and MAPLE 
routines. Then by using the algebraic manipulation capabilities of MAPLE, we searched 
for the point where d'^f^^^/dj'^ = at j = as numerically we know this is where the 
critical point is. The exact values are given in table El as rationals, together with an eight 
digit approximate decimal value. The simple values obtained reflect the simple nature of 
the model. By way of comparison, in the spin-one model these fractions of integers are 
replaced by ratios of polynomials in exp{— a} where a is the extra physical parameter of 
the spin-one model. The values of these integers for the Ising model results depend on the 
details of the lattice used. We have no explanation for the fact that the same numbers 
appear twice in the first column, once in the numerator, once in the denominator, though 
related behaviour is known for the LDE expansion of the partition function, at least in 
zero-dimensions (i.e. the simple integral of exp{— x^ — Ax^} for x from — cxo to -(-cxo). 

Note the absence of a ID case in table[T] This is because the method does not produce 
any critical points in the ID case, which is physically correct — there are no phase 
transitions in the ID Ising model. We have already alluded to quantitative similarities 
between our model and mean field theory. The fact that we do not find a phase transition 
in the ID case is a step better than mean field theory, which does predict a phase transition 
in the ID Ising model. 

Pinally, we also used MAPLE to check the critical behaviour in the Ising model. Using 
the expressions for three dimensions we were able to show that the critical exponent 6 is 
indeed precisely 3. 

32 



Order 


Kjq 


for 2D 




1 


1/4 






2 


4/12 = 


1/3 = 


.3333333333 


3 


12/(104/3) = 


9/26 = 


.3461538462 


4 


(104/3)/92 = 


26/69 = 


.3768115942 


5 


92/(3608/15) = 


345/902 = 


.3824833703 


6 


(3608/15)/(9148/15) = 


902/2287 = 


.3944031482 


7 


(9148/15)/(485788/315) = 


48027/121447 = 


.3954564543 




r\jc. 


for 3D 




1 


1/6 






2 


6/30 = 


1/5 = 


.2000000000 


3 


30/148 = 


15/74 = 


.2027027027 


4 


148/706 = 


74/353 = 


.2096317280 


5 


706/(16804/5) = 


1765/8402 = 


.2100690312 


6 


(16804/5)/15746 = 


8402/39365 = 


.2134383335 


7 


15746/(7742666/105) = 


826665/3871333 = 


.2135349762 




Kjq 


for 4D 




1 


1/8 






2 


8/56 = 


1/7 = 


.1428571429 


3 


56/(1168/3) = 


21/146 = 


.1438356164 


4 


(1168/3)/(7976/3) = 


146/997 = 


.1464393180 


5 


(7976/3)/(272176/15) = 


4985/34022 = 


.1465228382 


6 


(272176/15)/(614728/5) = 


34022/230523 = 


.1475861411 


7 


(614728/5)/(262409816/315) = 


4840983/32801227 = 


.1475854242 




K 


^ for d 




1 


l/{2d) 






2 


(2d)/(4rf2 _ 2d) = 


l/(2d-l) 




3 


(4^2 _ 2rf)/(8d3 - 8^2 + 4rf/3(?)) 







Table 5: Exact critical points for the 2D, 3D and 4D Ising models The values are obtained 
by using the linear variational parameter j only. For 3D the Monte Carlo result is 
fi;c = 0.221654 [13 IZ21 ESI • 

5 Conclusions 



We have studied the LDE together with several different optimisation schemes for the 
Ising, model, spin-one model and full A0^ field theory. In doing so we have gone four 
orders higher than any previous LDE lattice study except for the pure gauge study of 
Kerler and Metz |^. In any case, to the best of our knowledge, the Ising and spin-one 
models have not previously been studied on the lattice using LDE. To reach our high 
orders we had to automate the whole diagrammatic expansion procedure. While similar 
expansions (high-temperature etc.) have been done to much higher orders elsewhere 
(e.g. see |SZ1EH]), these non-LDE studies do not have the optimisation stage of the LDE 
method. This means the limitations in terms of computing are quite different for LDE, 
we must calculate the whole expression for the free energy many many more times and 
our expressions have several extra parameters. Thus we have achieved one of our aims 
which was to demonstrate that high orders of LDE lattice expansions can be calculated 
without using significant computing resources. In doing so we have indicated some of the 
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key issues and our solutions to several problems. 

It is clear to us that the method could be improved in several places and higher orders 
reached, e.g. through the use of the 'free diagrammatic expansion' [77] . Though we could 
never match the orders produced for straight high-temperature expansions, experience 
with LDE in QM has been that the variational aspect more than compensates for the 
additional work or in our case, low numbers of terms in our expansions. 

We studied several different optimisation schemes, PMS (Principle of Minimal Sensi- 
tivity) and FAC (Fastest Apparent Convergence) but in this particular case, only our use 
of the minimum free energy principle seemed to work. It is hard to find detailed com- 
parisons of optimisation schemes or to find the use of the minimal free energy condition 
in the literature of LDE on the lattice, so it is hoped that our results have moved the 
debate on. 

We also note that all our work, and that of the LDE lattice literature, is exclusively 
on a simple cubic lattice, though of various dimensions. It would be straightforward to 
repeat our calculations for other lattices, just as high temperature expansions have been 
done. We merely note that calculations up to third order on a hypertetrahedral lattice 
have shown the expected variation in the position of critical points but there is no sign 
of a change in the critical exponents. 

We set out to use the calculation of critical exponents as a test of the LDE method, 
at least for this lattice approximation. As such we have shown that the LDE method 
has failed to capture the non-perturbative infra-red physics of scalar filed models near 
the critical point as all our results have given the incorrect mean field values. This is 
despite us pushing the calculation to far higher orders and studying the region close to the 
transition extremely carefully, taking great care to keep numerical errors under control. 
The fact that low orders produce these mean field values is not too surprising, as the 
exact analysis of the Ising model demonstrates, but we are unable to explain why the 
LDE method does not start to move towards the correct values as we study higher orders. 
In most studies where exact results are known (typically quantum mechanics) moving to 
higher and higher orders in the LDE method rapidly brings one to good approximations 
of non-perturbative results. Our conclusion is that this form of LDE on the lattice is 
unable to access the non-perturbative aspects of field theories near the critical point. 

This is not necessarily the end of the use of LDE for field theories. We are studying 
a theory which is not-asymptotically free, using a space-time lattice and only looking at 
one type of physics. LDE calculations for QCD, using continuous space-time methods^'^ 
or even for other quantities might show better results. One might try to find a different 
way of applying the LDE in this case. After all, the high-temperature expansion gives a 
series which itself necessarily has no transitions yet by using Fade methods and similar 
one can extract remarkably good results. Such methods are not directly applicable to 
LDE as it does not calculate coefficients of a series. LDE only provides a sequence of 
results which in some testable cases converge extremely fast to the correct answer. For 
the LDE the most obvious variation is to try different ansatz for the trial action 5*0. We 
have indicated that it must be ultra-local for calculational reasons, but beyond that we 
have great flexibility For instance we tried ansatz such as 

Lo=j(P + k<p' + q\<P\^ + l<P^ (82) 

with non-polynomial power p a variational parameter. However our studies at low orders 



"'^^We note, however, that that the results of Gandhi and McKane [3D] for critical exponents using 
LDE with continuum space-time methods were only slightly better. See also the work of Ogure and Sato 
[SUES], andChiku [21]. 
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showed no promise and we did not pursue this further. 
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A Growing unique configurations 

For diagrams of n links (an abstract graph), we could easily construct all possible con- 
figurations (embeddings of this graph on a space-time lattice) by throwing down n links 
in all possible ways on all the lattice graph edges in an (n + 1)*^ size box, but every con- 
figuration then appears in {n + 1)*^ copies differing only by translation invariance. This 
is clearly highly inefficient at the orders and dimensions we work at so the aim of the 
algorithm outlined here is to produce a complete set of configurations of a certain size, 
once and only once, modulo translation symmetry. The algorithm is inspired by those 
described for percolation simulations j66|l67j. 

The starting point is to define some relation on the edges in the configurations, so that 
we can always decide that any certain line is 'greater' or 'less' than any other line. To do 
this, we first have to define a relation on the points of the lattice, and so we say that a 
point X = {xi, X2, ■ ■ ■ , Xd) is greater than a point y, i.e. x > y, ii the first coordinate of 
X is greater than the first coordinate of y, xi > yi- If the two coordinates happen to be 
equal, then we compare the second coordinates in the same way. If these turn out to be 
equal, we continue moving on to the next coordinate until we have exhausted all the d 
coordinates that label a point on the lattice. 

X > y <^ 3 i s.t. Xi > Hi and Xj = yj ^ j < i, i,j G {1,2, ... ,d} (83) 

Of course, if all the coordinates turn out to be the same, the two points pi and p2 are 
equal. 

In our configurations edges, have no directions but it is always useful to list them 
with the least vertex first. An edge is then an ordered pair of two (neighbouring) points 
e = {x,x^a),x < ajnn- To define a relation between two edges e and f, we say that e 
is greater than f if e's first vertex is greater than the first vertex of f . If the two points 
happen to be equal, then we compare the second coordinate in the same way. 

e> f ^3i s.t. Ci > fj and Cj = fj \/ j < i, i,j G {1, 2} (84) 

The relation we defined on the lines will ensure that every configuration has a unique 
least line. This link goes through the least vertex of the configuration which we choose 
to be the origin O on the space-time lattice. All the configurations we create have the 
origin as the least vertex as then we are guaranteed that all our configurations are not 
related by translation symmetry to any of the others we produce. 

What then remains to be done is to ensure that we grow all possible configurations 
once and only once from this lowest vertex origin. A recursive algorithm which does 
this is as follows. For an existing configuration of n links, C„, the set of possible lines 
to be added is an ordered list of links, £„, in which each link satisfies the following five 
conditions: 

1. they never contains a vertex less than the origin O (eliminates translations), 
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2. they start from an existing vertex in the configuration (we are not interested in 
disconnected graphs), 

3. they do not he on top of any other hnks (we are not interested in non-simple graphs), 

4. they are always greater than the least link in the existing configuration C„ (ensures 
we visit each configuration only once, and none are related by translation invariance 
as the first and least line remains the least line) 

5. each link has not been investigated at previous orders even if it is not in the existing 
configuration C„ (ensures each configuration is visited once and only once) 

We then add a link to C„, starting with the least link in the list of possible links £„, 
creating a configuration C„+i of ra + 1 links. We then recurse, considering possible links 
to add to this n + 1 link configuration, Cn+i- We stop recursing when R links have 
been added. Once we have considered adding one link to the n-link configuration, say 
Cj G Cn, and all the possible larger configurations coming from that, before we consider 
the next possible link, Cj+i G £„ for this n-link configuration C„, we mark the link just 
investigated, Cj, in the lattice as having been used. This enables the implementation of 
the fifth condition in the production of the lists Cm {m > n) and is vital to eliminate 
repetitions of the same configuration. 

Note we do not generate configurations in the order of their links. The relation ()84|1 
gives a natural order for the links in every configuration but many configurations have 
links which are connected to each other only via links which are greater than both of 
them. This is why the fourth condition only requires potential new links to be greater 
than the first one. 

To see how this works, or to test the voracity of alternative algorithms, the example 
oi d = 2 space-time dimensions on a square lattice up to configurations of i? = 3 links is 
often sufficient. For instance consider the two size 2 configurations of figure IT31 With the 
first coordinate up the page, the second to the right of the origin O, we exclude the lattice 
vertices shown as open circles by rule one and so eliminate repetitions related by space- 
time translation. It also follows that configuration (a) is generated before configuration 
(b) so this is why the link labelled X has can not be used in the latter, by rule five. 
Without this limitation, the order 3 configuration using link 4 in case (a) would have 
been reproduced by case (b) with link X since the latter link is legitimate by all the other 
rules. 



B Free Energy Density 

The calculation of the free energy density is based on the formalism developed in section 
12] with the free energy density given as an expansion of order R 

/(^) = -EV- = -EV'^ (85) 

n=0 n=0 

We also know that /o = Ib-Zq, and for the other /„'s, we have to make use of all the 
different connected graphs of size n and their multiplicities. For n = 1 this is a particularly 
easy task, but the n = 2 term is illustrative, yet simple enough to work through as an 
example. 
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Figure 15: Two different second order configurations are shown in solid lines. The first 
coordinate is measured up from the origin (9, the second coordinate is measured to the 
right of the origin. Thus the vertex points with open circles are never used to ensure 
that the origin O is always the smallest vertex, thus eliminating space-time translations 
of diagrams. The links which can be added to produce third order diagrams are given as 
dashed lines, numbering reflecting the edge ordering as defined in the text, i.e. Cj < Cj if 
i < j. Edge X in case (b) is crossed out to indicate it is excluded by our rule five ensuring 
that when link 4 is added to configuration (a), the resulting order 3 configuration is not 
produced from configuration (b). 



The diagram counting is the first step in the numerical evaluation of any /„. The C/C++ 
program implementing the algorithm for counting multigraphs without loops is used to 
find the multiplicities of all the graphs of size 2 and less. Note that the multiplicities 
will depend on the dimension of the problem, i.e. the dimension of the lattice that the 
configurations are generated on. For definiteness we will consider the example of a d = 3 
cubic lattice. The program will output the results in a text file, and for our example we 
get 

3 • — . 3 #=» 30 • — . — • (86) 



The numbers preceding the graphs are the multiplicities, the number of distinct ways 
these graphs can be embedded on a space-time lattice up to translations. Thus /ii,/22 
are 
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The cumulant expectation values then have to be expanded in terms of statistical averages 
using equation (jHHj) . For example. 
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'>c = (' 



■>■-(— X("X 



(89) 



These graphs can be in terms of fields, physical and variational parameters, and the 
expectation values with respect to the Li trial Lagrangian fields at one site. Using (jHUj) 
these last contributions are expressed in terms of the Jp integrals of ()31|). For instance 



\ _ Ji Ji 



(90) 
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However, the J„ integrals contain 6i in their exponents through Li. Terms such as fu 
are aheady order (^2)^ so they can have their J„ integrals simply replaced by /„ of (fT^ . 
However fu have only one power of delta so they require us to expand the J„ integrals 
to first order in 61 and the terms from /lo = In(Jo) require full second order expansions. 
For the Ising model case we have 

J„ := /„ + <52AJ/„+i + ^(52AJ)2/„+2 + 0{{62)') (91) 

where AJ = J—j. Inserting into (jHTj) truncating at order of delta R = 2 and then setting 
the (5's equal to one, we finally get the following expression for second order term of the 
full free energy, /2: 

h = 3o[(0»^>,-(0>:]+3[(0^>^(0>; 

+12 [<0>^ [U - J)(0^>^ + ik- a){<p\ + (/ - g){<p\] 

- Wo [0- - ^) Wo + (^ - «) Wo + (^ - ^) Wo] 

+2 [{j -J){k- a){<j>')^ + (j -J){1- g){<l>')^ + (k-a){l- ^)(0^> J 

-0- - ^f Wo - (^ - « W>o - (^ - ^W>o 

-2[(j -J){k- «)(0>,(0^>„ + (j -J){1- gMM")o 

+ (fc-a)(/-^)(0^>^(0^>J (92) 

where {<jf)^ = I J I,. 

We have shown the above expression in its entirety to demonstrate the general form 
of the /„. Thus we see that the calculation of /„ is broken down into a sum of terms 
(the number of which will get rather large as we go to higher orders) depending on the 
physical parameters J, a and g, the variational parameters j, k and /, and also the field 
expectation values (0^)q- Note that the (0^)„ actually depend only on the variational 
parameters. From a numerical perspective, however, we look upon the (0^)^, as quantities 
in their own right since their evaluation is a non-trivial task. Thus, mathematically we 
write f^^\j,a,g;j,k,l) to indicate the variables which contribute to the calculation of 
y(-R)_ On the other hand, numerically we recognise that in order to calculate f^^\ we also 
need to address separately the calculation of a set of (0^) 's. 

Maple is the ideal tool to perform all the necessary expansions, starting from the 
appropriate list of diagrams and multiplicities such as (J86|l . through to obtaining the 
expressions for the /„. For the particular case of n = 2, this will be an expression of the 
form shown in equation ()92|1 above. Maple will output /2 as code in the C programming 
language, optimising the sum for minimum computation. We note that this output 
was always checked up to 3rd order against code that we had available from the study 
of a complex (p^ theory [IH]. By performing this 3rd order check routinely, we gained 
confidence in the validity of our diagram counting and code generating techniques. 

C Numerical Accuracy 

The standard double precision available to C/C++ (type double) is limited to around 17 
decimal places of precision on floating point numbers. ^^ The nature of the calculations 

^"^Technically, this will depend on the computer architecture and the particular C/CH — \- implementa- 
tion. We used a standard PC based on an AMD Athlon XP 2100+ processor (a 32-bit architecture), and 
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in our model is such that there are many additions and subtractions, with a very wide 
range of values. Here we have to deal with the common problem of terms in a summation 
conspiring to cancel out in such a way that smaller terms fail to contribute to a sum in 
which their contribution is, in fact, important. This is the rounding problem. 

C.l Free Energy Density 

Now consider an example from our actual model — the calculation of f^^^ for a set 
of physical and variational parameters. The numerical values of the parameters in this 
example are taken from one of the points visited by the program as it searches for the 
critical point. This is a realistic scenario for the double precision program, since it is 
taken at a stage where the program is exploring variations of the order of 10"^^ in the free 
energy. This is the sort of accuracy we would expect to achieve from an implementation 
using C/C++'s standard double precision. The values of interest are listed in table IHl 



fn 


Smallest 


Largest 


Total 


n= 1 
n = 2 
n = 3 


1.4997840002 ■ 10-^° 
4.6741039631 • 10^22 
2.0210522783 • lO^^s 


2.7472168610-10-1° 
7.4509483795 ■ 10"^ 
7.4509483795 ■ 10"^ 


-9.1453475295-10-11 
-7.4509483756 ■ 10"^ 
-6.4570642135 • lO^n 



Table 6: The table summarises the absolute values of the smallest and largest terms in 
the sums which constitute the /„ terms, for all orders up to the third. The total value of 
each order is given too. These particular values are taken from a real calculation of f^^^ 
close to the critical point. 

The third order proves already to be numerically both interesting and challenging. 
The smallest term being of the order of 10~^^ and the largest of 10~^. The total result 
is of the order of 10~ii, and is arrived at by summing 93 termsi^ whose values will range 
from the smallest to the largest quoted in the table. It is obvious that if we are limited 
to double precision floating point arithmetic, we could easily end up completely losing 
contributions from terms that might actually be important, especially, one suspects, near 
physical transition points. One could argue that since the largest term is of the order of 
10"^ and the total is about 10~ii, we are actually only in danger of losing approximately 
9 decimal places in this case, which should still be well covered by the 17 decimal places 
that double precision offers. 

However it is not the free energy itself but differences in the free energies for similar 
parameter values which we are interested in, as we are trying to find a set of variational 
parameters which minimise the free energy. For instance, the values discussed above are 
derived from a point which is used to calculate differences in free energies of size 10~i^ 
only. Requiring greater accuracy of the free energy would require calculations involving 
even smaller values of the variational parameters, which also means a greater range of 
values appearing in the terms that contribute to the /„. 

Even if this level of accuracy is acceptable for this 3rd order case that we are present- 
ing, higher orders will further increase the disparity of the terms. For all these reasons 
we decided to use an arbitrary precision package, the GNU Multiple Precision (GMP) 
[7T] library. GMP fixes a number of bits which are allocated to the mantissa of a floating 
point number. This means that by setting the precision of the mantissa to, for example. 



we used Linux with the GNU Compiler Collection v3.2 implementation of C/C- 
^^In the 4D case. 
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128 bits, the largest number we can accurately represent is approximately 2^^^ ^ 3.4-10^®. 
In other words, we will have approximately 38 decimal places of accuracy. When we dis- 
cuss the results we will find numerous occasions where we really use the extra precision 
made available to us by the GMP library. For instance, we were able to minimise the 
free energy density to one part in 10"^''. 

C.2 Integration 

It has been emphasised that the calculation of the free energy density (and for that 
matter, any other quantity with the LDE formalism) will depend in particular on field 
expectation values at a single site taken with respect to the trial action So which depends 
on the variational parameters j, k and /, and of the form 

So{j,k,l)=j<P + k<j)' + l<p^ (93) 

We introduced in (fTH|l the notation Ip for the general form of the relevant integrals 

/ + 00 
d0 (jf exp [-J0 - A;02 _ /04j ^94^) 

-oo 

These are trivial for the Ising and spin-one models but non-trivial for the full (p^ model. 
This is the case we consider in this section. 

Being quartic in the exponential, the calculation of the /p's is a necessary numerical 
consideration. Numerical integration comes with its own set of problems though. For 
small j, the most interesting case for transitions, and for p odd, the integrand is nearly 
odd so to avoid severe rounding problems we used the half-range integral forms 

r+oo 

Ip = 2 d(f)(jfcosh{j(l))e~'">''-^*' for p even (95) 

Jo 

r+oo 

Ip = 2 d00Psinh(-j0)e-'='^'-''^' for p odd (96) 

Jo 

Having reduced the /p's to integrals over half the real line, we have (approximately) 
halved the number of interpolations that need to be performed in order to numerically 
evaluate the integrals. This will be an obvious computing time bonus, but we have also 
given ourselves a break from the dangerous 'rounding-off' error that appears in sums 
when terms are of opposite signs. 

Another important issue is the infinite range of integration. Our integrand is sup- 
pressed at large (p by the negative terms in the exponential. We will typically have 
I = 25, and, for example, if we were to introduce the cutoff for the upper limit at = 3, 
the exponential would look like exp{— 3j — 9k — 2025}. Thus the variational parameters 
j and k would have to be significantly large negative numbers to create a dent in the 
assumption that this exponential is vanishingly small, and in practice they are of the 
same order of magnitude as the corresponding physical parameters and so 0(10) at most. 
So it is easy to find a suitable cutoff value for in the Ip integrals. The numerical tests 
discussed below make it clear that we can indeed trust our numerical integration. 

Note that for every single evaluation of the free energy in the 0^ model, which will 
depend on the physical parameters as well as the variational parameters, we will have 
to evaluate a complete set of /p's, from p = through to whatever the upper limit on 
p turns out to be. For example, in a 3rd order expansion we are looking at the range 
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< p < 12, the upper limit coming from a term in which the quartic 0^ from the 
exponential is cubed upon expansion. For a 7th order expansion this range will thus go 
up to 28. We used a standard Romberg integration method, with additional polynomial 
extrapolation of successive refinements to zero stepsize jTO] . However we adapted this to 
perform the integration of all the required Jp's in parallel. Thus, instead of working with 
a single integrand (corresponding to a specific value of p), we manipulate a whole vector 
of elements, one entry for the integrand at each value of p required. 

The integrals and the numerical methods used to calculate them are not unusual. 
However, we did require much higher accuracy than normal so we will note how we 
tested our integration. For the purpose of testing, we used the following values for the 
variational parameters: ^^ 

j = 8.16343354579642367932991357923 ■ 10~^° (97a) 

k = -7.6 - 1.58180882338368170019383733935 • 10"^ (97b) 

/ = 25 (97c) 

Using 256 bit floating point arithmetic, approximately 77 decimal places of precision, 
table [7| shows a comparison of our results for the above integrals for various values of 
p with the best results we could get for the same integral calculated using numerical 
integration in Mathematica |78j . 



p 



h 



1.67769083077124890028663394770701193593 
1.67769083077124890028663394770701193592521428271904249117805060996740 



1.81907694824780353539804027376 ■ 10-^° 
1.81907694824780353539804027376452421473264747834415743104674406673865-10-1° 



0.0143833783157673970638326143289412430245 
0.01438337831576739706383261432894124302449287542710830784439069879080 



11 



5.606129663619230401281915955874438188097 • lO^^^ 
5.60612966361923040128191595587443818807152405548720830361149104778581- lO^^^ 



12 



0.00068673673058880728698436218620528130020 
0.00068673673058880728698436218620528130019933560098958463918835374638 



Table 7: The table shows results of evaluating the integrals Ip for various p, with the 
variational parameters as given in ()97p . For each value of p we show the result from 
Mathematica (first row), and the result from our integration routine (second row). Our 
results are calculated with 256 bit floating point arithmetic; we show approximately 70 
decimal places (we expect the results to be that accurate). 

A slightly different check was performed on the Iq integral for j = k = and / = 1 
and / = 25, which is simply 2r(5/4)(Z)~i/^. Mathematica has no problem evaluating the 
Gamma function to very high precision. We found that our integration routine matched 
the 150 decimal places of Mathematica's results. i''' 

^^The apparent randomness of the values is due to the fact that these constituted a point of contention 
in the early stages of development of the program; it subsequently assumed a status of a testing point. 
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To achieve this accuracy we used 512 bit floating point arithmetic. 
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Arithmetic / Tolerance 


Kjq 


64 bit / 10-1° 


0.2027041103 


64 bit / 10-1^ 


0.202702704077629 


128 bit / 10-1° 


0.2027041103 


128 bit / 10-15 


0.202702704077255 


128 bit / 10-20 


0.20270270272418008017 


128 bit / 10-25 


0.2027027027027865987084316 


128 bit / 10-30 


0.202702702702702784632396012135 


256 bit / 10-20 


0.20270270272418008017 


256 bit / 10-25 


0.2027027027027865987084316 


256 bit / 10-30 


0.202702702702702784632395797442 


256 bit / 10-35 


0.20270270270270270302274056635402966 


256 bit / 10-^0 


0.2027027027027027027039528506075906985909 



Table 8: The table summarises the critical inverse temperatures k^. found for the 3D Ising 
model at 3rd order. Each row represents a run with a different precision on the internal 
arithmetic and tolerance value for the simplex optimisation. The values of Kc are given to 
as many decimal places as is allowed by the tolerance required on the minimum of f^^\ 

D Benefits of High Precision Arithmetic 

Many studies based on LDE compatible methods (e.g. ISOllSIllSnilZn]) noted the similarity 
between their critical exponents and those obtained by mean field theory. Indeed, in 
preparing the results for |45j, we also noticed that critical exponents calculated with our 
3rd order LDE method produce values equivalent to mean field theory results. For this 
reason, one of the main goals of this study was to estimate critical exponents with high 
numerical accuracy, in particular to establish any possible deviations from the mean field 
results as higher orders of the LDE are considered. 

As a prerequisite, we must find the critical point. Table |H1 summarises the estimates 
of Kc for the 3D Ising model at 3rd order. The values were found by searching for the 
value of K where {^0) (k, J = 0) switches between zero and non-zero values. We will 
now explain the significance of the results given in the table. 

The results are all obtained with an implementation using the GNU multiple precision 
library (GMP). The first two rows show results obtained using 64 bits of precision, which is 
comparable to the built-in double precision arithmetic of C/C++. The remaining rows were 
calculated with significantly higher precision. 128 bit arithmetic gives us approximately 
38 digits of precision, while for 256 bit arithmetic we get approximately 77 digits. Every 
row is also labelled by a certain tolerance, which is the termination criterion for the 
simplex minimisation routine. 

To understand what the numbers are telling us, we have to explain the route by which 
they are obtained. The quantity which is minimised to the tolerance setting is the third- 
order free energy density /*^3) of the free energy. The /p's are evaluated to high accuracy 
since they are simple functions in the Ising model. Once evaluated, these are combined 
with the physical and variational parameters to form Z^^-*. The minimisation routine 
keeps the physical parameters {k and J) constant, but varies the variational parameters 
(only j in this case) until it hits a minimum in f^^\ The procedure will continue minimi- 
sation until it is satisfied that it has reached the tolerance which is required of it. That is 
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why we quote k^ in table |H1 with as many decimal places as the free energy was minimised 
to. 

Although /^^^ is minimised to within some quoted tolerance of the minimum point, 
there is no reason to assume that the variational parameters have been fixed to the same 
accuracy. In fact, experience shows that relatively large variations in the variational 
parameters (of the order of e, say) produce relatively small variations in /^^^ (of the 
order of e^). This should come as no surprise, because we expect linear variations not 
to contribute at a stationary point. Moreover, this lends support to the reasoning of the 
PMS criterion, since it is our goal to fix the variational parameters in a way which will 
leave the physics as independent of the variational parameters as possible. 

For this simple Ising model case the exact location of the critical point is 15/74 k, 
0.2027027027. . . (see section 1^?^ . We see in table |H1 that the numerical accuracy in Kc 
increases as we increase the tolerance as long as this is paired with enough accuracy 
provided by the internal arithmetic. If we require a tolerance of 10~^° on the free energy, 
we are confident that j is correct to within approximately 10~^° of its value. This explains 
why the values of k^ in table |H1 follow the exact result up to about half of the decimal 
places displayed. 
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